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Abstract 


This  report  constitutes  an  integrated  geophysical  investigation  of  the  lithospheric 
structure  in  the  region  of  Libya.  It  is  separated  into  five  sections.  The  final  four  sections 
represent  papers  that  will  be  submitted  to  different  scientific  journals  for  publication.  In 
the  first  part  of  the  study,  we  utilized  a  seamless  mosaicking  approach  based  on  the 
commercial  Environment  for  Visualizing  Images  (ENVI)  software  package  to  create 
mosaics  of  two  geologically  interesting  portions  of  Libya.  In  this  study  we  presen  :  a  step 
by  step  method  of  mosaicking  Landsat  4  satellite  images.  Firstly,  we  performed 
histogram  matching  to  give  images  the  same  color  scale,  then  we  used  a  cutline 
feathering  technique  to  blend  suture  areas  and  finally  we  overlaid  the  images  to  form  the 
two  mosaics.  The  resulting  mosaics  were  then  combined  with  structural  features  and  the 
seismicity  map  of  the  area.  The  resulting  mosaics  were  proven  to  be  useful  in  identifying 
recently  active  faults  and  shows  great  potential  for  verification  of  other  faults  and  in 
natural  hazard  assessment. 

For  the  second  portion  of  my  research,  we  made  use  of  over  6,000  Free  Air 
corrected  gravity  data  in  conjunction  with  other  geological  and  geophysical  data  to 
develop  a  3D  density  model  for  northern  Libya.  We  used  a  gravity  modeling  program 
(SURFGRAV)  to  develop  the  3D  density  model  by  manipulating  it  to  accurately  predict 
large  areas  of  Free  Air  anomaly  shown  in  the  data.  The  residual  gravity  anomaly  values 
were  calculated  by  subtracting  predicted  Free  Air  anomaly  from  the  observed  Free  Air 
anomaly.  The  results  were  satisfactory  for  uplifted  areas  of  Libya  while  there  were 
significant  mismatches  in  basin  areas.  The  density  model  was  iterated  and  used  as  a 
starting  model  for  the  final  portion  of  the  study. 

We  then  used  the  Nafe-Drake  relationship  along  with  other  geological  data  to 
convert  the  3D  density  model  to  a  3D  velocity  model  (LIBYA3D)  for  the  region.  Two 
earthquakes  having  source  receiver  paths  sampling  much  of  the  modeled  area  were  used 
to  perform  ID  and  1.5D  validation  tests,  and  the  results  were  compared  to  those  from 
previous  studies.  The  results  also  showed  that  the  new  3D  velocity  model  is  valid  and 
superior  to  the  global  model.  However,  until  there  is  sufficient  earthquake  data  acquired, 
and  we  are  able  to  perform  2D  and  3D  modeling  we  may  not  be  able  to  see  the  true 
improvement  of  LIBYA3D  as  compared  to  the  other  regional  models. 

In  the  final  part  of  this  report,  we  utilize  a  3-D  geophysical  model,  “LIBYA3D,” 
in  the  context  of  regional  discrimination  to  investigate  small  to  moderate  size  natural 
events  (m  3.7  -  5.7)  in  the  Libyan  region  occurring  from  May  1990  to  May  2000. 

b 

Observations  of  Pn  and  Sn  energy  at  varying  frequency  ranges  (0.5  -  8  Hz)  and 
propagation  paths  recorded  by  broadband  instrumentation  led  us  to  investigate  the 
attenuation  affects  of  the  geologic  structure  of  the  Libyan  region  and  the  adjacent 
Mediterranean  Sea.  We  perform  multi-frequency  analysis  to  discern  which  frequencies 
have  the  optimum  signal  to  noise  ratio  (SNR)  for  these  regional  phases.  Pn  and  Sn  SNR 
show  stronger  amplitudes  at  higher  frequencies  (4-8  Hz)  along  propagation  paths  where 
continental  crust  is  limited.  Conversely,  we  only  observe  strong  Pn  and  Sn  amplitudes  at 
lower  frequencies  (0.5  -  4  Hz)  along  propagation  paths  where  continental  crust  is 
prevalent.  We  utilize  LIBYA3D  to  develop  relationships  in  our  data  as  a  function  of 
crustal  path,  Moho  variance,  and  volume. 
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Section  1 


Executive  Summary 

Our  final  report  is  in  four  parts.  Part  one  is  the  text  of  a  paper  to  be  submitted  for 
publication  entitled  “A  Seamless  Mosaicking  of  Landsat  IV  Images  of  Northern  Libya 
and  its  Application.”  In  this  paper  we  utilized  a  seamless  mosaicking  approach  based  on 
the  commercial  Environment  for  Visualizing  Images  (ENVI)  software  package  to  create 
mosaics  of  two  geologically  interesting  portions  of  Libya.  The  resulting  mosaics  were 
proven  to  be  useful  in  identifying  recently  active  faults  and  show  great  potential  for 
verification  of  other  faults  and  in  natural  hazard  assessment. 

Part  two  is  the  text  of  a  paper  to  be  submitted  for  publication,  entitled  “The  Use  of 
Free  Air  Gravity  Data  to  Construct  A  3D  Density  Model  for  the  Region  of  Libya.”  In  this 
paper  we  made  use  of  over  6,000  Free  Air  corrected  gravity  data  in  conjunction  with 
other  geological  and  geophysical  data  to  develop  a  3D  density  model  for  northern  Libya. 
The  results  were  satisfactory  for  uplifted  areas  of  Libya  while  there  were  significant 
mismatches  in  basin  areas.  The  density  model  was  iterated  and  used  as  a  starting  model 
for  the  next  two  portions  of  the  study. 

Part  three  is  the  text  of  a  paper  to  be  submitted  for  publication,  entitled 
“Development  and  Validation  of  a  3D  Velocity  Model  for  the  Libya  Region.”  In  this 
paper  we  used  the  Nafe-Drake  relationship  along  with  other  geological  data  to  convert 
the  3D  density  model  to  a  3D  velocity  model  (LIBYA3D)  for  the  region.  We  use  two 
earthquakes  having  source  receiver  paths  sampling  much  of  the  modeled  area  to  perform 
ID  and  1.5D  validation  tests.  The  results  showed  that  the  new  3D  velocity  model  is  valid 
and  superior  to  the  global  model. 

Part  four  is  the  text  of  a  paper  submitted  for  publication  in  the  Bulletin  of  the 
Seismological  Society  of  America  entitled  “Analysis  of  High  Frequency  Propagation  in 
the  Libyan  Region  Utilizing  A  3D  Geophysical  Model”.  In  this  paper  we  utilize  our  3D 
geophysical  model,  “LIBYA3D,”  in  the  context  of  regional  discrimination  to  investigate 
small  to  moderate  size  natural  events  (m  3.7  -  5.7)  in  the  Libyan  region  occurring  from 

b 

May  1 990  to  May  2000.  Observations  of  Pn  and  Sn  energy  at  varying  frequency  ranges 
(0.5  -  8  Hz)  and  propagation  paths  recorded  by  broadband  instrumentation  led  us  to 
investigate  the  attenuation  effects  of  the  geologic  structure  of  the  Libyan  region  and  the 
adjacent  Mediterranean  Sea.  These  empirical  observations  could  be  used  to  aid  future 
explosion  discrimination  work  in  the  Libyan  region  if  the  need  arises. 
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Section  2 


A  SEAMLESS  MOSAICKING  OF  LANDSAT  IV  IMAGES  OF  NORTHERN 
LIBYA  AND  ITS  APPLICATIONS 


Summary 

The  general  goal  of  mosaicking  is  to  create  a  single  seamless  image  by  aligning  a 
series  of  overlapping  images.  The  result  is  a  larger  image  that  portrays  a  large  area  of 
interest.  In  this  paper,  we  utilized  a  seamless  mosaicking  approach  based  on  the 
commercial  Environment  for  Visualizing  Images  (ENVI)  software  package  to  create 
mosaics  of  portions  of  northeastern  and  northwestern  Libya.  The  mosaics  were  created 
from  thirteen  georeferenced  Landsat  4  satellite  images.  The  two  areas  chosen  (NE  Libya 
and  NW  Libya)  are  geologically  interesting  in  terms  of  seismic  activity,  complex  surface 
geology,  and  several  styles  of  active  faulting.  The  resulting  mosaics  have  great  potential 
to  aid  in  the  verification  of  previously  mapped  faults,  identification  of  other  unmapped 
faults,  and  assessment  of  natural  hazards  and  risk. 


Introduction 

Mosaicking  is  a  technique  of  image  processing  that  was  initially  used  to  manually 
combine  aerial  photographs.  Satellite  image  mosaicking  is  the  process  of  digitally 
combining  several  individual  satellite  images  so  that  the  boundaries  between  the  original 
images  are  not  seen.  Image  mosaicking  requires  geometric  corrections,  so  that  all  of  the 
images  match  in  a  geographic  sense.  When  transformed  to  any  local  coordinate  system 
(e.g.  UTM)  and  then  mosaicked,  the  result  is  an  orthophoto  map  sheet  (Afek  and  Brand, 
1998).  The  transformation  of  the  images  to  a  local  coordinate  system  (Moik,  1980)  is  also 
called  orthorectification.  For  this  study,  the  images  to  be  mosaicked  were  orthorecti fied 
to  UTM  zone  34. 

Image  mosaics  can  be  used  for  different  applications  such  as,  construction  of 
virtual  environments  and  virtual  travel  (Kim  et  al.  2003,  Szeliski  1996,  and  Szeliski  and 
Shum  1997);  for  scene  stabilization  (Hansen  et  al.  1994,  and  Irani  et  al.  1995a,  b);  and  for 
scene  change  detection,  video  compression  and  video  indexing  (Irani  and  Anandan 
1998).  Here  we  use  it  for  the  purpose  of  fault  identification,  verification  and  natural 
hazard  mapping. 

Creating  an  image  mosaic  was  originally  performed  using  photomechanical 
devices  (Mullen,  1980;  Afek,  1998).  Creating  a  mosaic  in  this  manner  has  always  been 
tedious  and  incurs  considerable  room  for  error.  Non  digital  mosaicking  is  still  in  use  (e.g., 
Vickers,  1993).  In  this  paper  we  propose  a  step  by  step  digital  method  of  mosaicking. 
This  method  is  not  fully  automated  but  it  is  a  quick  and  relatively  simple  technique  which 
will  provide  a  way  of  creating  mosaicking,  especially  in  cases  where  a  remote  sensing 
software  is  not  equipped  with  a  built  in  automated  mosaicking  routine. 

Mosaicking  remotely  sensed  images  for  the  northern  region  of  Libya  will  assist  in 
the  analysis  of  structural  features  in  the  area.  In  this  paper  we  utilize  a  seamless 
mosaicking  technology  to  create  mosaics  of  northern  Libya.  First,  we  performed 
histogram  matching  based  on  the  portions  of  adjacent  scenes  that  overlapped.  We  then 
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used  a  cutline  feathering  technique  to  blend  the  suture  areas.  Finally,  we  overlaid  the 
images  to  form  a  complete  mosaic. 

Mosaicking  satellite  images  will  provide  insight  into  recent  faulting  and 
seismicity  in  the  region.  This  paper  represents  part  of  an  effort  to  better  understand  the 
ways  in  which  seismic  waves  propagate  across  the  Libyan  lithosphere  so  that  earthquakes 
can  be  discriminated  from  nuclear  testing.  The  results  of  this  project  will  be  useful  to  the 
nuclear  monitoring  community  and  anyone  interested  in  locating  seismic  events  within 
Libya  or  in  modeling  the  propagation  of  seismic  waves  across  North  Africa. 

Study  Areas  and  Data 

Northern  Libya,  situated  in  the  Mediterranean  foreland  of  the  African  shield, 
extends  over  a  platform  of  cratonic  basins  and  includes  classic  examples  of  continental 
rifts  and  domal  uplifts.  These  basins  were  active  as  recently  as  the  Late  Tertiary  and 
Holocene  (Goudarzi,  1980).  Faults  with  various  trends  are  present  in  Libya,  but  the  major 
fault  systems  trend  sub-parallel  to  the  Red  Sea.  This  study  is  part  of  a  large  effort  to 
better  understand  the  lithospheric  structure  beneath  Libya,  utilizing  the  seismic  studies  as 
well  as  satellite  imagery  to  identify  recently  active  faults. 


Figure  1.  Basemap  of  the  study  area  showing  the  area  covered  by  the  remotely 
sensed  Landsat  IV  images.  The  labels  indicate  the  path/row  of  the  images.  The  red 
lines  indicate  faults  cutting  rocks  of  Mesozoic  age  and  younger,  while  the  black  lines 
indicate  faults  cutting  rocks  of  Precambrian  age. 
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Thirteen  remotely  sensed  Landsat  4  images  covering  two  different  regions  of 
Northern  Libya  were  carefully  selected  for  this  study  ().  Four  of  the  images  are  182/38 
(path/row),  183/38,  184/38  and  183/37,  covering  a  portion  of  northeastern  Libya 
(Cyrenaica  Uplift).  This  area  is  the  second  most  seismically  active  portion  of  the  country 
and  is  the  site  of  much  recent  seismic  activity.  The  geology  of  this  area  is  mostly 
composed  of  Triassic  rocks  and  Quaternary  surficial  deposits.  A  few  major  basement 
high  angle  reverse  faults  are  located  in  the  area.  The  other  nine  images  ( 1 86/39,  1 86/40, 
186/41,  187/38,  187/39,  187/40,  187/41,  188/38,  and  188/39)  cover  the  northwestern 
portion  of  the  country.  This  region  is  the  most  seismically  active  portion  of  the  country, 
highly  faulted,  geologically  much  more  complex,  and  contains  a  portion  of  the  Hun 
Graben.  The  Hun  graben  trends  NW-SE  and  is  a  classical  example  of  well-exposed 
normal  block  faulting  (Suleiman  and  Doser,  1995).  In  the  northwest  area,  strike-slip 
faulting  is  common  with  a  general  northwesterly  trend.  Analysis  of  these  images  will  help 
in  identifying  the  more  recently  active  faults  and  previously  unmapped  faults. 

In  order  to  best  interpret  these  faults,  it  is  necessary  to  mosaic  the  individual 
images.  The  mosaicked  image  will  also  serve  as  a  base  map  for  overlaying  other 
information  such  as  earthquake  locations,  which  will  help  interpreters  to  better 
understand  the  relationship  between  earthquakes  and  faults. 

Seamless  Mosaicking  Technology 

Mosaicking  is  the  art  of  combining  multiple  images  into  a  single  composite 
image.  It  can  be  used  to  combine  pixel-based  images,  or  as  a  means  for  combining 
georeferenced  images  into  a  single  image  covering  a  larger  geographic  area.  The 
seamless  mosaicking  approach  based  on  the  commercial  Environment  for  Visualizing 
Images  (ENVI)  software  package  used  in  this  study  was  first  developed  by  Xie  (2002) 
and  improved  in  this  study.  This  approach  involves  three  procedures:  (histogram 
matching,  cutline  feathering,  and  ordering  overlain  images)  that  are  discussed  in  turn 
below. 

Histogram  Matching  Based  on  Overlay  Parts 

In  order  to  match  two  images  that  will  be  mosaicked,  Stepl  is  to  mark  the 
overlapping  parts  of  the  two  images  as  a  region  of  interest  (ROI),  for  example,  the  images 
with  marked  ROIs  are  shown  in  Figure  2.  Step  2  (Figure  3)  is  to  match  the  two  images  by 
matching  the  histograms  of  the  two  ROIs.  The  histogram  of  the  image  184/38  is  set  to 
stretch  type  arbitrary,  since  its  output  will  be  modified  to  resemble  the  histogram  of 
image  183/38.  The  output  histogram  of  the  base  image  (183/38)  is  then  dragged  and 
dropped  onto  the  output  histogram  of  image  184/38  and  applied.  The  procedure  is 
repeated  twice  for  the  other  colors  of  the  RGB  color  spectra.  This  gives  the  overlapping 
region  almost  exactly  the  same  grayscale  distribution  (Figure  4).  The  improvement  of  this 
step  compared  to  other  histogram  matching  is  that,  instead  of  matching  the  entire  images, 
we  only  match  the  overlapping  areas  (ROIs).  There  is  a  risk  when  matching  the  entire 
images  due  to  large  black  areas  (DN  is  0)  outside  the  image.  These  large  zero  value  areas 
significantly  impact  the  distribution  of  Display  Number  (DN)  value,  and  then  impact  the 
histogram  matching.  Figure  4  is  the  result  of  histogram  matching. 
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Feathering  Technique  to  Blend  Suture  Areas 

After  finishing  the  histogram  matching,  a  cutline  feathering  technique  was  used 
for  blending  the  suture  areas.  A  cutline  must  be  defined  using  the  annotation  tool  prior  to 
mosaicking,  usually  along  the  seam  where  one  desires  the  mosaicked  boundary  to  be 
matched.  The  annotation  file  must  contain  a  polyline  (white  line,  image  184/38  of )  as  a 
cutline  that  is  drawn  from  edge-to-edge  and  a  symbol  (green  star)  that  is  placed  in  the 
region  of  the  image  that  will  be  cut  off  ().  The  cutline  distance  is  used  to  create  a  linear 
ramp  that  averages  the  two  images  across  that  distance  from  the  cutline  outwards. 


Figure  2.  The  overlapping  parts  of  the  two  images  are  marked  as  ROIs  (red  on  left 
image  and  green  on  the  right  image).  The  left  image  is  path/row  of  184/38  and  the 
right  image  is  183/38.  It  is  clear  that  the  right  image  is  significantly  brighter  than 
the  left  image. 
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Figure  3.  Illustration  of  the  histogram  matching  method.  Notes:  The  histogram  of 
the  ROI  in  183/38,  the  base  histogram,  is  used  to  match  the  histogram  of  ROl  in 
184/38.  The  histogram  for  image  183/38  remains  unchanged,  but  for  image  184/38, 
the  output  histogram  more  closely  resembles  the  output  histogram  in  image  183/38. 


Figure  4.  Images  after  histogram  matching. 


Ordering  overlain  Images 

It  is  important  to  place  the  images  in  a  particular  order  before  trying  to  decide  which 
image  will  have  the  cutline  and  symbol  on  it.  Here  are  some  basic  principles  that  need  to  be  kept 
in  mind: 

1)  The  base  image  will  always  stay  on  the  bottom,  and  the  other  adjacent  images  that  will 
be  cut  off  while  mosaicking  will  always  stay  on  the  top. 

2)  The  cutline  must  be  within  the  bounds  of  both  the  images  being  cut  and  the  base 
image,  otherwise  there  will  be  a  gap  left  in  the  final  mosaicked  image. 

3)  The  cutline  needs  to  be  assigned  for  every  two  adjacent  images,  otherwise  there  will  be 
a  seam  left  along  the  bounds  of  the  top  image  because  there  is  a  linear  ramp  area  created  along  a 
cutline  by  averaging  the  two  images  using  a  specified  cutline  distance  (300  pixels  in  our  case).  If 
a  cutline  is  not  defined  along  two  adjacent  images,  there  will  be  no  calculation  along  the  suture 
area  and  a  seam  will  be  left  along  the  bounds  of  the  top  image. 
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Results  and  Summary 

Figure  5  shows  the  final  mosaicked  image  from  4  individual  scenes  for  the  northeastern 
portion  of  Libya.  From  this  image,  we  cannot  see  any  visible  effects  along  the  seams  and  color 
contrast  from  image  to  image,  which  means  that  this  is  a  good  mosaic.  In  Figure  6  we  have 
overlain  mapped  Precambrian  faults,  Mesozoic  faults  and  seismicity  of  the  area  over  Figure  5. 
This  mosaic  clearly  indicates  the  correlation  between  faults  and  possible  changes  in  vegetation.  It 
is  noteworthy  however,  that  the  images  depict  surface  reflectance  that  is  not  always  related  to 
vegetation.  Having  the  faults  overlain  on  top  of  the  mosaicked  image  helps  to  discern  what  a 
fault  looks  like  on  a  satellite  image.  This  helps  us  in  identifying  suspect  or  unmapped  faults.  We 
may  also  be  able  to  determine  relative  age  of  some  of  the  faults  based  on  their  geomorphic 
expression  and  crosscutting  relationships  with  other  features.  Since  the  faulting  styles  exhibited 
in  the  Cyrenaica  Uplift  and  the  Hun  Graben  area  are  so  different,  we  may  even  be  able  to 
identify  the  styles  of  some  of  the  unmapped  faults  based  purely  on  the  satellite  images. 


Figure  5.  Final  mosaicked  image  for  the  northeastern  portion  of  Libya. 


Another  overlay  including  earthquake  locations  (Figure  6)  may  enable  us  to  associate 
earthquakes  with  the  specific  causative  faults. 

From  the  mosaicked  images  several  unmapped  and  probably  unknown  topographic 
features  have  been  noted,  particularly  on  the  mosaicked  image  in  the  northwest  region  of  Libya. 
The  northwest-southeast  trending  faults  that  define  the  flanks  of  the  Hun  Graben  are  clearly  seen, 
but  their  exact  position  may  be  slightly  off  due  to  pixel  shifting  (e.g.,  Scheik,  2004).  Additional 
active  faults  can  also  be  seen  in  the  area  based  on  the  similarities  shown  on  the  mosaicked  image 
to  other  mapped  faults  (Figure  7).  The  topography  of  the  Cyrenaica  area  has  been  outlined  and 
can  be  used  to  produce  more  accurate  geographical  maps  of  the  region,  therefore  improving  our 
ability  to  produce  geo-hazard  zonation  maps. 
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Figure  6.  This  flgure  shows  the  two  mosaicked  area  overlain  by  mapped  faults  in  Libya. 
The  blue  faults  are  faults  that  cut  Precambrian  and  younger  rocks  (Goudarzi  et  al.,  1978), 
the  red  faults  cut  only  Mesozoic  (Goudarzi  et  al.,  1978)  and  younger  rocks,  and  the 
magenta  faults  are  other  mapped  faults  of  unknown  age  in  the  region  (Hearn  et  al.,  2001). 
The  red  dots  show  the  seismicity  of  the  region  (1975-2004).  Earthquakes  are  of  magnitude 
3.5  and  above  (NEIC  PDE  catalog)  (see  Table  1). 
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Figure  7.  This  figure  shows  the  nine  mosaicked  images  for  northwestern  Libya.  Part  a) 
shows  the  image  after  mosaicking  and  part  b)  shows  the  image  overlain  by  seismicity, 
mapped  faults  and  assumed  unmapped  faults  (green  lines). 


Conclusions 

Figure  6  and  7  show  both  mosaics  overlain  by  mapped  faults  of  known  age  in  Libya 
along  with  recent  seismicity  of  the  region  (Table  1).  For  the  mosaicked  section  of  northeast 
Libya,  it  is  clear  that  all  seven  seismic  events  occurring  in  the  region  are  associated  with  the 
Precambrian  aged  faults  that  outline  the  borders  of  the  Cyrenaica  uplift.  Faults  cutting  rocks  of 
Mesozoic  age  and  younger  are  less  extensive  and  are  therefore  less  obvious  on  the  mosaicked 
image.  For  this  reason  we  may  not  be  able  to  distinguish  faults  of  Mesozoic  age  from  the 
mosaicked  satellite  image.  However,  on  zoomed  images  with  directional  filters  applied,  we  were 
able  to  identify  linear  features  that  align  with  mapped  faults  that  appear  to  be  surface  expressions 
of  other  unmapped  faults. 

The  results  were  similar  for  the  mosaicked  image  for  the  northwest  region.  The  sharp 
change  in  reflectance  in  the  center  of  the  image  demarcates  the  flanks  of  the  Hun  Graben,  along 
with  the  dominantly  northwest  southeast  trending  Precambrian  faults  in  the  region.  Focal 
mechanism  studies  by  Suleiman  et  al.  (1995)  suggest  there  are  high  angle  faults  in  this  region. 
Again  the  faults  of  Mesozoic  age  and  younger  were  not  extensive  enough  to  be  used  to  predict 
the  location  of  other  faults  of  the  same  age  on  a  macroscopic  scale.  The  sharp  change  in 
reflectance  in  the  southern  portion  of  the  image  shows  the  presence  of  a  mapped  fault  of 
unknown  age.  The  principle  of  cross  cutting  relationships  could  not  be  applied  to  discern  the  age 
of  this  fault  as  it  was  quite  isolated  from  other  features  of  known  age.  The  seismicity  of  the 
northwest  portion  of  Libya  seems  to  be  mostly  associated  with  faults  of  Precambrian  age.  Three 
out  of  four  of  the  seismic  events  in  the  northwestern  region  are  located  along  a  mapped 
Precambrian  fault,  suggesting  that  these  may  be  near  vertical  faults  if  they  are  linked  to  the 
earthquakes.  The  single  event  that  is  located  to  the  east  of  the  nearest  mapped  fault  would  also 
suggest  that  the  fault  is  dipping  to  the  east  and  they  are  linked  at  depth. 
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Table  1.  List  of  events  occurring  in  the  study  region  (1975  -  2000).  Only  events  between 
latitude  18N  -  33N  and  Longitude  25E  -  10E  are  listed  in  this  table.  (Source  NEIC  PDE 
catalogue). 


Date(yy  mm  dd) 

Origin 

Time(sec.) 

Lat.  (deg.) 

Long. (deg.) 

Depth 

(km) 

Mag.  (mb) 

75 

01  01 

07:03:51.40 

32.46 

21.21 

39 

4.2 

75 

09  26 

22:42:21.40 

28.65 

17.24 

33 

78 

06  07 

07:21:21.50 

25.80 

15.09 

33 

4.1 

83 

10  17 

22:39:07.67 

31.29 

15.99 

10 

4.4 

87 

06  28 

00:50:17.86 

32.82 

24.36 

24 

5.4 

87 

09  13 

23:55:00.01 

32.64 

24.37 

33 

4.1 

88 

01  28 

15:48:08.15 

32.41 

21.15 

10 

4.8 

88 

01  28 

19:12:17.11 

32.24 

21.16 

10 

4.4 

90 

05  18 

18:27:51.27 

31.71 

24.80 

10 

4.3 

92 

07  01 

15:13:32.38 

31.47 

15.76 

10 

4.3 

93 

09  09 

00:32:40.17 

32.37 

20.65 

33 

4.5 

96 

10  30 

10:43:02.17 

32.33 

20.50 

10 

4.4 

99 

02  10 
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Section  3 


THE  USE  OF  FREE  AIR  GRAVITY  DATA  TO  CONSTRUCT  A  3D  DENSITY  MODEL 
FOR  THE  REGION  OF  LIBYA 


Summary 

We  present  a  new  3D  density  model  for  the  Libya  region  that  includes  the  topographic 
relief  and  five  individual  subsurface  layers,  the  deepest  being  at  a  mantle  depth  of  50  km.  The 
model  was  derived  utilizing  over  6  000  gravity  measurements  in  conjunction  with  other 
geophysical  and  geological  data  obtained  from  oil  exploration,  published  research,  and  other 
academic  institutions.  To  develop  the  density  model,  we  used  “SURFGRAV”,  a  newly 
developed  general  purpose  gravity  modeling  program.  To  evaluate  the  validity  of  our  model,  we 
computed  the  predicted  Free  Air  gravity  anomalies  for  the  study  area  and  compared  our  results  to 
the  observed  Free  Air  gravity  anomalies.  The  gravity  modeling  showed  that  our  new  density 
model  is  able  to  predict  large  portions  of  the  observed  Free  Air  anomalies,  especially  in  uplifted 
areas.  With  the  exception  of  the  northwestern  and  northeastern  regions  where  data  are  sparse  and 
thinning  of  the  crust  would  improve  the  fit  the  residual  Free  Air  gravity  is  less  than  50  mGals, 
indicating  that  a  more  detailed  model  is  required  to  produce  a  better  match  between  observed  and 
calculated  gravity  values. 


Introduction 

The  tectonic  evolution  of  Libya  has  yielded  a  complex  crustal  structure,  which  is 
composed  of  a  series  of  basins  and  uplifts  that  began  forming  in  the  Lower  Paleozoic  (Goudarzi, 
1980).  Tectonic  deformation  has  continued  to  the  present  time,  with  continued  basin  subsidence 
and  eruption  of  lavas  associated  with  the  Tibesti  Mountains  in  southern  Libya.  A  considerable 
amount  of  oil  exploration  has  been  undertaken  in  the  area  and  numerous  studies  have  been 
published  on  the  shallow  (<10  km  depth)  geology  and  geophysics  of  the  region.  In  addition,  over 
6  000  gravity  measurements  are  available  for  the  northern  Libya  region. 

To  achieve  our  goal  of  creating  a  3D  density  model  of  the  study  area,  we  primarily  used 
SURFGRAV,  a  general  purpose  gravity  modeling  program  that  bases  its  computations  on 
representing  the  earth  by  gridded  Digital  Elevation  Maps  (DEMs)  with  an  associated  grid  with  a 
map  of  geologic  units  (Baker,  2001).  In  this  approach,  topography  and  near  surface  geology, 
greatly  simplified  in  conventional  gravity  interpretations,  are  included  in  the  model  to  determine 
detailed  structure  lying  below  the  surface.  Veilleux  (1999)  first  utilized  this  technique  in  a  study 
of  the  shallow  subsurface  of  the  Pajarito  Plateau  region  surrounding  Los  Alamos,  New  Mexico, 
which  is  an  area  of  extreme  topography  consisting  of  surficial  and  near  surface  low  density 
volcanic  tuffs.  Using  this  technique,  Veilleux  (1999)  was  able  to  better  delineate  the  eastern  and 
western  boundaries  of  the  Velarde  graben  located  beneath  the  Pajarito  Plateau. 

In  this  study,  we  use  available  gravity  data,  in  conjunction  with  other  geologic  and 
geophysical  control,  to  construct  a  3D  model  of  density/geology  for  northern  Libya  and 
surrounding  regions.  We  then  applied  a  DC  shift  to  our  calculated  Free  Air  results  in  order  to 
obtain  a  range  of  gravity  values  that  were  comparable  to  the  range  in  observed  gravity.  We  then 
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compared  our  resulting  calculated  gravity  values  to  determine  the  regions  where  our  density 
model  appears  to  match  the  observations. 


Geology 

Libya  is  located  on  the  Mediterranean  foreland  of  the  African  shield.  It  is  a  platformal  area 
containing  cratonic  basins  and  classic  examples  of  continental  rifts  and  domal  uplifts  (Figure  8). 
Precambrian  igneous  and  metamorphic  rocks  comprise  a  small  fraction  (less  than  1%)  of  the 
surface  of  the  geology  of  Libya  (Suleiman,  1985).  There  are  basement  exposures  of  mainly 
granitic  igneous  rocks  and  some  metamorphic  rocks  in  south  central  Libya  in  the  Tibesti  area  of 
the  Nubian  uplift  (Figure  8  &  Figure  9).  Vachette  (1964)  used  radiometric  dating  to  show  that 
various  granitic  rocks  in  the  Nubian  uplift  were  formed  500  to  600  m.y.  ago.  The  Ben  Ghnema 
batholith  in  the  center  of  the  Nubian  Uplift,  northwest  Tibesti  is  a  complex  of  plutons  and  is 
considered  to  be  similar  in  many  ways  to  the  southern  California  and  Sierra  Nevada  bacholiths 
(Pegram  et  al.,  1976;  Ghuma  and  Rogers,  1978).  Ghuma  and  Rogers  (1978)  concluded  :hat  the 
Ben  Ghema  batholith  was  formed  by  the  subduction  of  oceanic  crust  westward  under  a 
continental  margin  located  in  the  Tibesti  area  during  the  Pan  African  orogeny.  Precambrian  rocks 
are  also  found  further  north  associated  with  the  Fezzan  and  Thiemboka  uplifts. 

Paleozoic  sedimentary  rocks  are  widespread  over  southern  Libya  (Figure  9).  These  are 
mainly  marine  shales,  siltstones,  limestones  and  continental  sandstones.  There  are  no  mapped 
Paleozoic  rocks  outcropping  north  of  latitude  28N.  However,  they  are  reported  in  the  subsurface 
in  many  drill  holes  in  the  northern  basins  (Suleiman,  1976;  Hammuda,  1980). 

The  Mesozoic  era  of  the  region  is  characterized  by  the  deposition  of  hundreds  of  meters  of 
sediments  in  shallow  and  deep  troughs  along  the  Mediterranean  shore  (Goudarzi,  1970). 
Sandstones,  shales,  limestones,  conglomerates,  and  various  undifferentiated  rocks  of  Mesozoic 
age  can  be  found  throughout  Libya. 

Rocks  of  Tertiary  age  (mainly  carbonates)  are  concentrated  in  the  northern  part  of  the 
country.  Pleistocene  and  Holocene  deposits  are  also  widespread  forming  the  Libyan  desert,  the 
dunes  areas,  the  gravel  plains  and  the  Mediterranean  coastal  plains  (Goudarzi,  1970).  A  great 
majority  of  the  Sirte  Basin  contains  Tertiary  and  Quaternary  deposits  (Figure  8). 

Tertiary  and  Quaternary  olivine,  basalts  and  phonolites  cover  a  large  portion  of  the  central 
part  of  the  country.  Piccoli  (1970)  reported  basalts  to  have  an  average  age  of  40  m.y.  However, 
an  average  age  of  45  m.y.  was  determined  for  the  southern  part  of  the  country. 

It  is  only  within  the  last  twenty  five  years  that  a  coherent  picture  of  the  plate  tectonic  history 
of  Libya  has  emerged.  This  is  partially  due  to  the  fact  that  two  distinctively  different  aspects  of 
tectonics  are  involved.  In  broad  terms,  the  evolution  of  Gondwana 


13 


12H 


16E 


20E 


8E 


24E 


32N 


2XN 


24N 


20N 


Figure  8.  Map  of  Libya  showing  the  spatial  extent  of  basins  and  uplifts  occurring  in  Libya 
(Based  on  Hearn,  2001). 

and  Pangaea  controlled  the  development  of  Libyan  tectonics  during  the  Palaeozoic  and  Early 
Mesozoic,  whereas  the  evolution  of  the  Tethys  and  the  Mediterranean  Seas  were  the  primary 
controls  on  the  tectonics  of  Libya  during  the  late  Mesozoic  and  Cenozoic. 
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Explanation  of  Legend 


C  Carboniferous  rocks 

CO  Carboniferous  and  Upper  Devonian  rocks 

Cf-P  Carboniferous  and  Lower  Trtassic 

D  Devonian  rocks 

DS  Devonian  Silurian  rocks 

Du  Middle  and  Upper  Devonian  rocks 

H  Holocene  roots 

J  Jurassic  rocks 

JTr  Jurassic  rocks  (Liassk  age) 

Jl  Jurassic  and  Tr iassic  rocks  (undivided) 

K  Cretaceous  rocks 
KC  Cretaceous  rocks  ( Cenomanian  age) 

KJ  Cretaceous  Jurassk 
Kl  Lower  Cretaceous 
Mi  Lower  Miocene  rocks 
Mz  Miocene  rocks 
N  Nubian  rocks  { lower  cretaceous) 


O  Ordovician  rocks 

OCm  Ordovkian  Cambrian  rocks 

P  Pal  eocene  rocks 

Pg  Intrusive  ignoeus  rocks 

Q  Surficial  deposits 

Qe  Quaternary  eolian  deposits 

Qp  Pliocene  deposits 

QT  Quaternary  Tertiary  deposits 

Qv  Quaternary  Vokanks 

S  Silurian  rocks 

T  Tr  iassic  rocks 

Ti  Pliocene  rocks 

TK  Upper  Cretaceous  rocks 

Tr  Traissk  Rocks 

Tr-J  Triassic  Jurask  rocks 

TrP  Upper  Paleocene  rocks 

pC  Precambrian  rocks 


Figure  9.  Generalized  geologic  map  of  Libya  and  adjacent  areas. 


The  main  ev  idence  for  the  evolution  of  the  northern  margin  of  Gondwana  has  come  from  studies 
on  the  ancient  cratons  and  the  intracratonic  mobile  belts  exposed  in  Algeria,  Chad,  Niger,  and 
Sudan  (Hallet,  2002).  Given  the  diverse  origins  of  North  Africa,  it  is  not  surprising  that  even 
after  the  break-up  of  Pangaea  Africa  has  not  operated  as  a  single  plate.  The  intraplate  weaknesses 
within  the  African  continent  have  had  a  significant  influence  on  the  tectonic  development  of 
Libya  (Hallet,  2002). 
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Previous  work 

Suleiman  et  al.  (1991)  carried  out  a  gravity  study  of  the  Sirte  Basin  in  Libya  using  gravity 
data  archived  at  the  Libyan  National  Oil  Corporation  and  available  from  major  oil  companies. 
They  produced  a  simple  Bouguer  gravity  map  for  the  Sirte  basin  which  was  used  for 
interpretation.  They  noted  an  increase  in  regional  gravity  values  from  south  to  north  with  a  rapid 
increase  at  the  coastline,  typical  of  the  continent-ocean  margin  due  to  thinning  of  the  continental 
crust.  They  found  that  the  most  prominent  gravity  anomalies  correlated  well  with  the  structural 
highs  and  lows  of  Libyan  basins  and  uplifts.  They  also  used  computer  modeling  to  construct 
three  gravity  profiles  (one  of  the  profiles  traversed  the  Murzuk  Basin,  the  Fezzan  Uplift,  the  Sirt 
Basin  and  the  Cyrenaica  Basin  in  Libya)  and  found  that  their  earth  models  confirmed  the  block- 
faulting  nature  of  the  basins  and  uplifts  and  indicated  thinning  of  the  crust  beneath  the  rifted 
areas. 


Dial  (1998)  utilized  gravity  and  seismic  waveform  modeling  techniques  to  determine 
models  of  lithospheric  structure  across  northern  Africa.  He  constructed  density  models  using 
gravity  data  constrained  by  previous  geological  and  geophysical  studies.  He  constructed  three 
lithospheric  profiles  based  on  Bouguer  gravity  modeling.  One  of  his  three  lithospheric  profiles 
traversed  Libya  (Figure  10  and  Figure  11).  He  also  used  wavelength  filtering  techniques  to 
investigate  the  relative  depth  and  extent  of  the  structures  associated  with  the  major  anomalies. 
He  noted  that  the  resulting  earth  model  showed  slightly  greater  crustal  thickness  in  the  Atlas 
Mountains  (WNW  of  Libya)  than  those  of  previous  studies  if  a  low  density  mantle  region  is  not 
included.  The  results  also  suggested  that  most  of  the  northern  African  mantle  is  fairly 
homogenous  in  density  with  the  exception  of  regions  located  beneath  the  Atlas  Mountains  and  to 
a  limited  extent  beneath  the  Hoggar  Mountains. 
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Figure  10.  Map  of  North  Africa  showing  the  location  of  gravity  stations  and  profiles 
constructed  by  Dial  (1998). 
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Profile  3 


Figure  11.  Dial’s  1998  model  for  gravity  profile  1  (see  Figure  10  for  location  of  profile).  The 
model  extends  for  2620  km  from  southern  Algeria  to  northern  Egypt  (Dial,  1998). 
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“SURFGRAV” 

Carrying  out  gravimetric  terrain  corrections  is  an  important  and  tedious  data  processing 
step  for  work  in  rugged  topography.  To  deal  with  this  problem  the  classical  method  of  employing 
fan-shaped  prisms  was  developed  by  Hayford  and  Bowie  (1912)  and  improved  by  Hammer 
(1939).  Bott  (1959),  Kane  (1962),  Nagy  (1966),  and  Plouf  (1966,  1977)  further  developed 
computer-oriented  terrain  correction  algorithms  which  represent  topography  as  rectangular 
prisms.  The  common  “draw-back”  related  to  these  methods  was  that  they  obtained  gravity 
terrain-correction  values  by  computing  volume  integrals  in  which  fan  shaped  or  rectangular 
prisms  are  used  to  approximate  the  topographic  relief  discontinuously  (Zhou  et  al,  1990).  Zhou  et 
al  (1990)  approximated  the  land  surface  using  dipping  triangular  elements,  as  opposed  to  the 
conventional  use  of  prisms,  and  they  found  that  this  triangular  element  method  is  more  efficient, 
accurate,  and  flexible  than  methods  used  in  previous  studies. 

The  3D  gravity  modeling  used  in  this  study  (SURFGRAV)  is  based  on  representing  the 
Earth's  surface  by  gridded  Digital  Elevation  Maps  with  associated  grids  representing  tops  and 
bottoms  of  geologic  units  (Baker,  2001).  The  program  requires  three  inputs  to  be  used  in  the 
modeling  process  (Figure  12).  They  are:  a  gravity  data  file,  gridded  surficial  and/or  subsurface 
geology,  and  a  digital  elevation  model  (DEM).  The  output  of  the  program  is  the  computed 
gravity  due  to  the  model  constructed,  which  is  presented  as  a  matrix  or  spreadsheet  data  set  with 
each  row  representing  a  gravity  station,  and  the  columns  containing  unit-density  gravity 
contributions  for  each  of  the  geologic  units  in  the  area  of  computation. 
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Figure  12.  Flow  chart  for  processing  and  modeling  3D  gravity. 


Computation  Techniques  Used  in  SURFGRA  V 

Existing  topographic  DEMs  are  a  poor  representation  for  terrain  corrections,  as  the  grid 
spacing  is  far  too  coarse  for  near  station  corrections,  and  far  too  fine  for  far  station  corrections. 
To  speed  the  process  up  to  the  point  where  the  computation  can  be  made  with  a  reasonable  effort 
and  accuracy  in  extreme  areas  of  topography,  a  number  of  geologic  body  representations  and 
numerical  approximations  are  used  by  SURFGRA V. 

First,  the  typical  geologic  body  is  broken  into  upper  and  lower  surfaces.  This  makes 
vertical  sided  prisms  easy  to  set  up  and  compute.  Second  the  stations  are  grouped  in  arbitrary 
bins  and  the  distance  from  the  bin  to  the  DEM/Geology  element  determines  the  calculation 
technique  used  by  the  algorithm  (Baker  2001). 

The  first  approximation  permits  a  geologic  classification  to  be  associated  with  the  surface 
of  the  geologic  body.  The  gravity  computation  for  the  upper  surface  of  the  geologic  body 
proceeds  by  accumulating  the  grid  defined  prism  contributions  into  their  respective  geologic 
classes  at  a  unit  density  value.  The  presumed  prism  extends  to  infinity  with  depth.  The  lower 
surface  contribution  can  be  accumulated  in  a  similar  way.  These  values  can  be  subtracted  from 
the  upper  surface  kernels,  giving  a  vertical-prism  type  variable  density  computation  (Baker 
2001). 


For  cases  where  the  size  of  a  prism  is  small,  compared  to  the  distance  of  the  station  from 
the  prism,  an  infinite  vertical  prism  is  represented  by  an  infinite  vertical  line  element  Within 
SURFGRAV,  all  grid/geology  elements  lying  beyond  a  specific  radius  (120  km  in  our  case)  are 
computed  using  the  vertical  line  element  formula.  The  120  km  radius  reflects  the  fact  that  it  is 
reasonable  to  use  this  approximation  for  regions  lying  at  a  radius  that  is  more  than  100  times 
greater  than  our  grid  element  size  of  0.87  km. 

In  a  SURFGRAV  setup  file  one  label  specifies  a  list  of  radii  in  kilometers  at  which 
various  approximations  are  used.  The  radius  at  which  one  wishes  to  invoke  these  approximations 
strongly  depends  on  the  granularity  of  the  geology  at  the  scale  of  the  modeling.  In  this  case,  we 
set  the  point  mass  radii  at  120  km,  the  row  mass  at  100km  and  the  adjusted  vertical  line  source  at 
1  km.  Large  blocks  of  homogeneous  geology  with  high  topographic  relief  require  larger  radii 
than  small  blocks  of  inhomogeneous  geology  on  a  flat  surface  (Baker  2001). 

Gravity  Data 

Over  6000  reduced  gravity  data  points  were  compiled  and  used  as  locations  at  which  to 
predict  gravity  values  in  the  modeling  process.  These  gravity  points  are  ultimately  used  in  this 
study  to  construct  a  3D  density  distribution  throughout  the  Libya  region.  The  majority  of  these 
points  were  collected  from  the  worldwide  gravity  database  maintained  by  the  National  Imagery 
and  Mapping  Agency  (NIMA)  which  is  now  called  the  National  Geospatial-Intelligence  agency 
(NGA).  The  other  points  were  taken  from  the  gravity  database  at  the  University  of  Texas  at  El 
Paso  (UTEP).  Even  though  there  is  a  large  amount  of  gravity  data  available  for  Libya,  the 
distribution  is  quite  uneven.  Approximately,  85%  of  the  gravity  data  available  for  Libya  are 
associated  with  the  Sirte  Basin.  This  is  due  to  the  high  level  of  petroleum  exploration  occurring 
in  this  region  of  Libya.  There  were  no  gravity  data  available  for  regions  of  Libya  located  to  the 
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south  of  latitude  25N.  In  his  previous  study  Dial  (1998)  had  obtained  additional  gravity  data 
supplied  by  GETECH  (associated  with  Leeds  University).  These  data  were  available  to  us,  but 
were  gridded  and  thus  no  station  elevations  were  available  so  they  could  not  be  used.  The 
distribution  of  the  gravity  data  used  in  this  study  is  shown  in  Figure  13. 


Figure  13.  Map  of  Libya  showing  the  locations  of  the  gravity  readings  used  in  this  study. 


Construction  of  the  Density  Model 

Topography 

The  Digital  Elevation  Model  used  for  the  surface  of  the  gravity  model  was  extracted  from 
three  main  sources:  Our  UTEP  gravity  database  which  had  gravity  station  elevations;  GTOPO30, 
which  was  released  by  the  USGS  (EROS,  1996);  and  the  Global  GIS  Database:  Digital  Atlas  of 
Africa,  USGS  Digital  Data  Series  (DDS62B)  (Heam,  2001).  These  data  sets  were  merged  and 
gridded  at  a  30  second  interval  to  form  a  new  combined  topography/bathymetry  data  set.  Figure 
14  is  a  3D  grid  portrayal  of  the  topography  of  Libya. 

Layers 

In  an  ideal  case,  each  geologic  unit  of  our  model  would  be  represented  as  a  separate 
polygon  with  an  associated  density  value.  If  various  sediment  thicknesses  within  basins  were 
known  or  could  be  calculated,  then  they  would  be  treated  as  individual  units  separate  and  apart 
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from  neighboring  units.  However,  these  distinctions  could  not  be  made  in  Libya,  as  no  data  with 
the  appropriate  resolution  are  presently  available  to  carry  out  this  exercise.  The  Mediteiranean 
Sea  was  treated  as  a  layer  within  our  model  with  the  thickness  estimated  from  bathymetry.  The 
standard  density  of  sea  water,  l.lg/cc,  was  assigned  to  this  water  layer.  We  proceeded  to  group 
the  surficial  deposits,  the  Cenozoic,  and  the  Mesozoic  units  into  a  single  layer,  thus  assuming  a 
uniform  density  down  to  the  Mesozoic/Paleozoic  interface.  The  densities  assigned  to  this 
combined  surficial  layer  were  derived  by  analyzing  the  geological  composition  of  the  outcrops 
and  assigning  their  average  empirical  densities  based  on  percentage  composition. 


Figure  14.  Vertically  exaggerated  DEM  of  Libya.  The  minimum  curvature  technique  was 
used  in  the  gridding.  The  gridding  was  done  at  30  sec.  spacing.  Data  used  to  create  the 
DEM  was  obtained  from  the  GTOPO30  database.  Elevation  is  in  meters. 


The  empirical  densities  (Table  2)  were  obtained  from  either  Carmicheal  et  al  (1989)  or  Dial 
(1998).  In  an  effort  to  reduce  the  number  of  geology  polygons  present  in  Figure  9,  neighboring 
polygons  of  geologic  bodies  having  densities  within  0.025  g/cc  of  each  other  were  combined  to 
form  a  single  polygon.  Figure  15  shows  the  resulting  average  density  map  of  the  upper  layer  of 
our  model. 
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Table  2.  Summarized  density  assigned  to  the  generalized  geology  of  Libya.  Table  also 
shows  the  age  and  basic  composition  of  the  different  units. 


Label 

Density 

Age 

Basic  Composition 

Cam 

2.60 

Cambrian 

course  sandstone,  quartzitic  cement 

Car 

2.55 

Carboniferous 

continental  sandstone,  and  siltstone 

PC 

2.70 

Precambrian 

quartzite,  phyllite,  marble,  schist,  gneiss 

Kct 

2.45 

Upper  Cretaceous  rocks 

limestone,  marl,  chert,  dolomite,  shale 

Kn 

2.50 

Lower  Cretaceous  rocks 

Sandstone,  conglomerate,  claystone 

Tmm 

2.20 

Mid  Miocene 

detrital  limestone,  calcareous  limestone 

Tv 

2.72 

Tertiary  Volcanic 

Olivine,  basalt,  phonolite 

Q 

2.00 

Quaternary  deposits 

Dunes,  loess,  alluvium,  calcarenite 

The  base  of  the  Mesozoic  and  the  top  of  the  Precambrian  basement  were  digitized  from  a 
1:  2  000  000  structure  contour  map  of  the  Libyan  Arab  Republic  and  adjacent  areas  (Goudarzi 
and  Smith,  1978).  Both  sets  of  contours  were  drawn  at  intervals  of  250  m.  In  northwestern  and 
northeastern  Libya  north  of  latitude  30N,  marine  Jurassic  and  Triassic  rocks  are  present  and 
serve  as  reliable  controls  for  the  mapping  of  the  base  of  Mesozoic  rocks.  Precambrian  rocks 
outcrop  in  southern  and  central  Libya,  again  serving  as  control  for  mapping  the  top  of  the 
Precambrian  basement.  The  base  of  the  Mesozoic  layer  and  the  top  of  the  Precambrian  basement 
defined  the  boundaries  of  the  Paleozoic  layer.  The  Paleozoic  layer  was  assigned  a  uniform 
density  of  2.55g/cc.  The  upper  boundary  of  the  lower  crust  was  set  at  a  depth  of  19  km,  which  is 
consistent  with  the  results  of  previous  studies  (Dial,  1998).  The  density  of  the  Precambrian  layer 
was  set  at  2.7g/cc  as  the  results  of  previous  gravity  models  dictated.  The  density  assigned  to  the 
lower  crust  is  3.0g/cc  and  its  upper  boundary  was  set  at  35  km  depth  and  its  lower  boundary  was 
defined  by  the  digitized  surface  of  the  Moho  (Dial,  1998)  (Table  3  and  Figure  16). 

Yousef  (1986)  had  previously  derived  a  crustal  thickness  map  for  northern  Africa  from 
surface  wave  dispersion  studies  across  Africa.  The  shear  velocity  model  computed  for  a  path 
across  North  Africa  Platform  and  the  Atlas  Mountain  suggests  an  average  crustal  thickness  of  37 
km  (Yousef,  1986).  Dial  (1998)  improved  upon  this  existing  crustal  thickness  map  using 
constraints  from  other  studies,  additional  petroleum  exploration  data  unavailable  to  Yousef 
(1986)  and  analysis  of  Bouguer  gravity  data.  The  major  improvement  made  by  Dial  (1998)  was 
to  thin  the  crust  to  38  km  in  north  central  Africa  near  the  Tibesti  uplift.  For  the  purposes  of  this 
study  we  used  a  density  of  3.3g/cc  (Dial,  1998)  for  the  upper  mantle.  The  maximum  depth  of  our 
model  was  50  km,  a  value  selected  to  be  greater  than  the  thickest  portion  of  the  crust.  Finally,  we 
used  a  30  second  grid  spacing  to  carry  out  our  density  modeling. 
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Figure  15.  Simplified  density  map  of  the  surface  geology  of  Libya. 
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Figure  16.  Illustration  of  the  arrangement  of  the  layers  used  in  the  density  model.  The  wavy 
lines  represent  the  irregular  surfaces  that  were  obtained  from  maps  and  other  sources.  The 
upper  surface  of  the  lower  crust  is  set  at  19  km  based  on  previous  studies.  The  base  of  the 
model  is  at  50  km. 
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Table  3.  Layers  used  in  the  density  model,  along  with  their  average  thickness,  assigned 
density  and  the  source  from  which  the  densities  were  obtained. 


Layer 

Ave  Thickness 

Density 

(g/cc) 

Source 

Mediterranean 

0.045  km 

1.10 

standard 

Surficial  Layers 

1.01  km 

2.00  -  2.75 

Carmichael  (1989) 

Paleozoic 

Precambrian 

1.057  km 

2.55 

Dial  (1998) 

Basement 

17.362  km 

2.70 

Dial  (1998) 

Lower  Crust 

15.362  km 

3.00 

Dial  (1998) 

Mantle 

15.957  km 

3.30 

Dial  (1998) 

Results/Discussion 

Isopach  maps  provide  a  convenient  means  by  which  we  can  visually  check  the  results  of 
the  gravitational  contribution  from  each  layer.  The  isopach  map  of  the  uppermost  layer  is  shown 
in  Figure  17.  The  surface  layer  is  thicker  in  basins  and  thinner  in  uplifts,  thus  we  expect  the 
gravitational  contribution  to  be  greater  in  regions  with  thicker  layers.  Figure  18  shows  the 
calculated  gravitational  contribution  for  the  uppermost  layer,  and  we  find  that  in  most  cases  it 
does  resemble  the  isopach  map  for  the  same  layer.  This  indicates  that  the  general  lateral  density 
variations  in  surficial  material  (Figure  15)  were  not  large  enough  to  overshadow  the  effect 
caused  by  changes  in  the  thickness  of  this  layer.  The  Sirte  basin  area  generally  returned  relatively 
high  calculated  gravity  contributions  and  so  did  the  Hamra  basin.  The  Cyrenaica  uplift  returned  a 
lower  than  average  calculated  gravity  contribution. 

Figure  19  shows  the  result  of  the  summation  of  the  calculated  gravity  contributions  of  all 
the  layers  to  a  depth  of  50  km.  Since  the  mantle  was  the  thickest  and  most  dense  layer  used  in  the 
study  (Table  3)  it  produced  a  gravitational  contribution  that  overshadowed  the  contribution  from 
the  overlying  layers.  In  order  to  meaningfully  compare  the  calculated  theoretical  gravity  with 
the  actual  Free  Air  anomaly,  we  applied  a  DC  shift  to  the  calculated  gravity  contributions.  The 
shift  was  selected  to  make  the  range  of  maximum  and  minimum  values  in  theoretical  gravity 
comparable  to  the  range  of  observed  data.  This  enables  us  to  subtract  the  observed  Free  Air  grid 
from  the  resulting  calculated  Free  Air  grid  to  obtain  the  residual  Free  Air  gravity  anomaly.  The 
resulting  anomaly  map  indicates  where  geologic  interpretations  may  have  been  misestimated. 
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Figure  17.  Isopach  map  of  the  upper  layer  of  the  model.  This  layer  combines  the  Cenozoic 
and  Mesozoic  layers.  Basins  appear  as  domes  since  the  layer  is  thicker  in  areas  where  the 
basins  are  Tilled. 


In  the  areas  of  mismatch  the  density  of  bodies  may  need  to  be  altered  or  the  thickness  of  bodies 
adjusted.  For  example,  in  areas  where  there  are  positive  residuals  it  implies  that  the  density  value 
assigned  to  a  geologic  unit  was  overestimated  or  the  thickness  of  the  unit  was  overestimated; 
conversely  negative  residuals  imply  underestimation  of  the  density  of  the  geologic  unit  or  the 
unit  thickness. 
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Figure  18.  Calculated  gravitational  contribution  of  the  uppermost  layer  of  the  model.  Data 
are  unavailable  south  of  latitude  26N,  hence  the  grid  is  truncated.  Note  the  similarities  in 
where  the  highs  and  low's  occur  as  compared  with  those  of  the  isopach  map  (Figure  17). 


29 


Figure  19.  Sum  of  the  gravitational  contribution  of  the  layers  used  in  the  model.  The 
contribution  from  the  mantle  had  the  most  significant  impact  on  the  results,  as  it  is  the 
thickest  and  most  dense  layer. 


The  residual  Free  Air  gravity  map  is  shown  in  Figure  20.  The  maximum  residual  was  296 
mGal  and  the  minimum  was  -150  mGal,  with  75%  of  the  residual  values  occurring  between  0  to 
1 36  mGal.  It  is  noteworthy  that  the  area  with  the  highest  density  of  data  (the  Sirte  Basin  area)  is 
not  necessarily  the  area  where  we  had  Free  Air  residuals  closest  to  0  mGals.  In  fact,  some  of  the 
largest  positive  residuals  occurred  in  the  Sirte  Basin.  This  could  be  attributed  to  the  fact  we 
assigned  a  single  density  to  the  entire  Sirte  Basin  region  down  to  the  top  of  the  Paleozoic  layer. 
The  basin  fill  may  not  be  as  consolidated  as  we  estimated  which  would  lead  to  us  overestimating 
the  gravity  values.  Also  our  modeling  may  not  have  adequately  portrayed  changes  in  the 
thickness  of  the  basin  fill,  especially  for  the  upper,  less  dense  portions  of  the  basin.  We  may  have 
also  overestimated  the  theoretical  gravity  values  in  the  Hun  Graben  and  the  Hamra  Basin  for  the 
same  reason.  The  areas  that  returned  the  more  suitable  density  correlation  seemed  to  be  those 
areas  that  were  uplifted,  (e.g.  Fezzan  Uplift  and  the  Cyrenaica  Uplift),  possibly  due  to  the  fact 
that  the  density  of  uplifted  areas  are  more  uniform  at  shallow  depth. 
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Figure  20.  Map  showing  the  residual  Free  Air  gravity  distribution.  The  residual  gravity 
values  were  calculated  by  subtracting  the  observed  Free  Air  gravity  from  the  calculated 
Free  Air  values.  Positive  residual  values  represent  areas  where  the  Free  Air  gravity  was 
over  estimated. 


Conclusions 

The  results  of  this  study  show  that  “SURFGRAV”  is  a  very  efficient  and  effective  tool 
that  can  be  used  to  estimate  gravitational  contributions  for  geologic  bodies.  We  accurately 
predicted  Free  Air  gravity  anomalies  for  the  uplifted  regions  in  Libya,  while  we  generally 
misestimated  anomalies  in  the  basins  in  the  neighborhood  of  1 50  mGals.  This  suggests  that  for 
the  basins  in  Libya,  we  need  additional  data  to  better  constrain  our  model.  With  more  data  we 
would  be  able  to  build  a  more  accurate  density  model  and  thickness.  It  is  still  possible  to  further 
refine  our  model,  but  at  present  it  is  the  best  density  model  of  crustal  Libya  available,  and  it  will 
serve  as  a  framework  for  deriving  a  3D  velocity  model  for  Libya. 
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Section  4 


DEVELOPMENT  AND  VALIDATION  OF  A  3D  VELOCITY  MODEL  FOR  THE  LIBYA 
REGION 

Summary 

A  detailed  velocity  model  for  the  Libya  region  will  play  a  vital  role  in  the  location  of 
earthquakes  and  other  events  within  Libya.  In  this  study,  we  developed  a  3D  velocity  model  for 
the  crust  and  upper  mantle  in  the  Libya  region  (called  LIBYA3D)  and  evaluated  its  suitability  for 
use  as  an  initial  model  in  tomographic  studies.  The  3D  velocity  model  of  Libya  was  developed 
by  integrating  the  results  of  previous  crustal  geophysical  studies  related  to  crustal  studies  in  the 
area  with  the  Nafe-Drake  relationship.  We  selected  two  seismic  events  occurring  in  the  region 
that  had  suitable  source  receiver  paths  across  the  study  area  and  performed  ID  and  1.5D 
validation  tests.  The  3D  model  is  accurate  to  within  1 .8  seconds  in  predicting  the  P-wave  arrival 
times  for  the  events  tested  and  was  also  useful  in  producing  regional  synthetic  waveforms  that 
more  closely  matched  the  actual  regional  data  than  other  velocity  models  for  the  region.  The 
results  of  tests  performed  on  the  3D  velocity  model  prove  that  LIBYA3D  will  be  a  particularly 
good  staring  model  for  future  studies  in  the  area,  which  may  include  wave  propagation,  regional 
discrimination  and  tomographic  studies. 

Introduction 

Confidence  in  the  ability  to  monitor  small  and  moderate  sized  events  is  crucial  for  the 
Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT)  (e.g.,  Gitterman  et  al.,  2003).  For  this  reason, 
it  is  important  to  have  a  3D  velocity  model  for  the  Libya  region.  Libya,  as  a  whole,  is  not 
considered  a  highly  active  seismic  region  (Gutenberg  and  Richter,  1965;  Campbell,  1968). 
However,  a  few  large  earthquakes  have  occurred  causing  large  scale  or  total  damage  (Kebeasy, 
1980).  Some  of  these  earthquakes  were  described  by  Minami  (1963)  and  Campbell  (1968).  3D 
velocity  models  are  usually  developed  using  tomographic  methods  that  require  the  initial 
velocity  model  to  be  a  reasonable  approximation  of  the  true  earth  structure  (e.g.  Johnson  and 
Vincent,  2002).  This  requires  extensive  data  and  station  coverage.  We  instead  derive  a  3D 
velocity  model  (hereafter  referred  to  as  L1BYA3D)  from  a  priori  geophysical  data  and  the  Nafe- 
Drake  (1957)  relationship  which  relates  density  of  geologic  bodies  to  velocity. 

LIBYA3D  was  derived  from  a  previous  study  conducted  by  Brown  et  al.  (2004)  in  which 
they  use  gravity  data  to  successfully  develop  a  3D  density  model  for  the  Libya  region.  Their 
density  model  consisted  of  a  layer  representing  the  Mediterranean  Sea,  the  Cenozoic  and 
Mesozoic  combined,  the  Paleozoic,  the  Precambrian,  the  lower  crust  and  the  upper  mantle 
(Figure  21).  VanDeMark  et  al.  (2004),  used  this  velocity  model  in  a  regional  discrimination 
study,  and  in  this  study,  we  use  the  Nafe-Drake  relationship  along  with  other  geophysical  studies 
conducted  in  the  area  to  convert  the  3D  density  model  to  a  3D  velocity  model  (see  Section  5). 
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Figure  21.  Map  showing  the  residual  Free  Air  gravity  distribution.  The  residual  gravity 
values  were  calculated  by  subtracting  the  observed  Free  Air  gravity  from  the  calculated 
Free  Air  values.  Positive  residual  values  represent  areas  where  the  Free  Air  gravity  was 
overestimated  (from  Brown,  2004). 

We  validated  our  model  by  conducting  several  tests.  First,  we  selected  seismic  events  that 
were  of  a  large  enough  magnitude  for  waveform  studies  to  be  performed.  The  ideal  events  would 
preferably  be  recorded  at  MEDNET  stations  and  have  suitable  source  receiver  paths  traversing 
our  area  of  study.  The  only  two  events  that  satisfied  this  criteria  were  event  921012  (yymmdd) 
(5.6mw)  and  event  980528  magnitude  5.5mw  (Figure  22).  In  this  paper,  we  present  the  results  of 
the  development  and  testing  of  our  3D  velocity  model  (LIBYA3D)  using  these  waveforms. 
Where  possible,  we  also  attempt  to  compare  our  results  with  other  proposed  velocity  models  for 
the  region. 
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Figure  22.  Seismicity  map  of  the  Libya  region  showing  earthquakes  as  red  diamonds, 
stations  as  blue  triangles  and  the  multicolored  lines  represent  source  receiver  paths.  Blue 
lines  show  source  receiver  pathways  that  were  suitable  for  testing  our  model.  Focal 
mechanisms  indicate  the  primary  events  used  to  validate  LIBYA3D. 


Existing  Geophysical  Work 

Intense  oil  exploration  activity  in  Libya  during  the  last  30  years  has  contributed  to  the 
acquisition  of  extensive  knowledge  of  the  geological  framework  of  Libya  and  its  continental 
shelf.  Few  regions  in  the  world  have  such  a  large  amount  of  geological,  geophysical  and  drilling 
data  available  that  have  been  collected  in  a  relatively  short  period  of  investigation. 

Kebeasy  (1980)  used  data  gathered  for  earthquakes  which  occurred  in  and  around  Libya 
for  the  period  262  A.D.  to  1977  to  investigate  the  nature  of  seismicity  and  seismotectonics  of 
Libya.  He  also  studied  the  distribution  of  earthquakes  in  time  and  space  and  calculated  the 
frequency  of  occurrence  of  shallow  earthquakes  in  the  past  40  years.  He  found  that  this 
distribution  of  earthquakes  in  time  and  space  suggests  classification  of  Libya  into  3  major 
seismic  zones  that  agree  well  with  the  distribution  of  major  tectonic  features.  He  also  found  that 
the  damage  due  to  earthquakes  covered  an  abnormally  large  area  around  the  epicenters  and 
attributed  this  mainly  to  the  small  attenuation  of  seismic  waves  generally  found  in  North  Africa. 

VanDeMark  et  al.  (2004)  have  utilized  LIBYA3D  in  the  context  of  regional 
discrimination  to  investigate  small  to  moderate  size  natural  events  (mb  3.7  -  5.7)  in  the  Libyan 
region  occurring  from  May  1990  to  May  2000.  Their  observations  of  Pn  and  Sn  energy  at 
varying  frequency  ranges  (0.5  -  8  Hz)  and  propagation  paths  recorded  by  broadband 
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instrumentation  led  them  to  investigate  the  attenuation  affects  of  the  geologic  structure  of  the 
Libyan  region  and  the  adjacent  Mediterranean  Sea.  Using  LIBYA3D,  they  performed  multi¬ 
frequency  analysis  to  discern  which  frequencies  have  the  optimum  signal  to  noise  ratio  (SNR) 
for  these  regional  phases.  Pn  and  Sn  showed  stronger  amplitudes  at  higher  frequencies  (4-8 
Hz)  along  propagation  paths  where  continental  crust  is  limited.  Conversely,  they  only  observed 
strong  Pn  and  Sn  amplitudes  at  lower  frequencies  (0.5  -  4  Hz)  along  propagation  paths  where 
continental  crust  is  prevalent.  They  further  utilized  LIBYA3D  to  study  relationships  in  their  data 
as  a  function  of  crustal  path  and  volume.  These  empirical  observations  could  be  used  to  aid 
future  clandestine  explosion  discrimination  work  in  the  Libyan  region  if  the  need  arises. 

Suleiman  and  Doser  (1995)  investigated  the  seismicity  of  Libya  in  order  to  better 
understand  the  earthquake  hazards,  the  tectonics  of  the  region,  and  the  present  day  evolution  of 
the  Sirt  Basin.  They  relocated  all  the  events  occurring  between  1931  and  1988  (Figure  23), 
determined  focal  mechanisms  from  first  motion  data,  and  used  waveform  modeling  to  determine 
the  source  parameters  for  the  1935  Mw  7  Hun  Graben  mainshock  and  its  two  largest  aftershocks. 
Their  study  showed  that  the  Hun  Graben  is  very  seismically  active.  They  showed  that  the 
majority  of  the  earthquakes  of  the  1935  sequence  occurred  along  the  eastern  side  of  the  Hun 
Graben  and  seem  to  align  very  well  in  a  NW-SE  trending  fashion  along  the  margin  of  the  graben. 
The  1939  sequence  was  found  to  be  located  in  the  Gulf  of  Sirt.  The  regions  to  the  NE  and  NW  of 
the  Sirte  Basin  are  seismically  active  (Figure  23)  with  almost  no  seismic  activity  occurring  in  the 
center  of  the  basin  itself. 

They  also  used  waveform  modeling  techniques  to  determine  source  parameters  for  the 
1935  main  shock  and  its  two  largest  aftershocks.  The  waveform  modeling  results  suggest  strike- 
slip  faulting  with  a  fault  plane  similar  to  the  strike  of  the  eastern  edge  of  the  Hiin  Graben.  The 
focal  mechanism  obtained  from  waveform  modeling  of  the  1935  mainshock  differed  from  that  of 
the  first  motion  studies  but  was  more  consistent  with  the  strike  of  faults  within  the  region.  A 
focal  depth  of  21  ±  3km  was  obtained  for  this  event.  More  recent  activity  (1939-1972)  appears  to 
have  extended  NE  from  the  graben  into  the  Mediterranean  Sea.  Suleiman  and  Doser  (1995) 
believe  that  if  the  zone  of  active  faults  associated  with  the  1935  and  1941  earthquakes  extends 
towards  the  north-west,  it  may  pose  a  significant  hazard  to  large  cities  in  the  region. 

Dial’s  (1998)  modeling  of  regional  and  teleseismic  waveforms  showed  that  seismic  wave 
propagation  throughout  the  Mediterranean  and  North  African  areas  is  complicated,  as  would  be 
expected  in  a  zone  of  plate  collision.  Dial  (1998)  made  a  comparison  between  the  velocity 
models  derived  from  his  study  and  previous  models  derived  in  Africa  and  the  Mediterranean  and 
found  that  the  upper  100  km  of  all  models  are  very  similar.  Observations  from  his  waveform 
modeling  led  him  to  conclude  that:  1)  either  a  pronounced  high-velocity  layer  lies  just  below  the 
asthenosphere  and/or  attenuation  of  the  Sn  phase  exists  throughout  much  of  the  Mediterranean 
and  northern  Africa  in  order  to  match  observed  waveforms  or  2)  the  1-D  reflectivity  modeling 
technique  used  in  his  study  (Randall  et.  al.,  1995)  was  insufficient  to  model  the  actual 
complicated  2-D  to  3-D  structure  of  the  area. 
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A  Earthquake  Location  (ISS  and  ISC  prior  to  1963  ) 
A  Earthquake  Relocation 

V  Faults  taken  from  Harding,  1984. 


Figure  23.  Epicenters  of  Libyan  earthquakes  (1931  -  1994).  Relocated  earthquakes  are 
denoted  by  solid  triangles  (from  Suleiman  and  Doser,  1995). 


Pasyanos  and  Walter  (2002)  attempted  to  estimate  the  crust  and  upper-mantle  seismic 
velocity  structure  in  North  Africa,  southern  Europe,  and  the  Middle  East  using  surface-wave 
dispersion  tomography  results  from  a  previous  study  (Pasyanos  et.  al.,  2001).  Their  surface  wave 
tomography  study  provided  high-resolution  coverage  across  the  region  for  more  than  6800 
Rayleigh  and  Love  wave  paths  over  waveform  periods  from  10-60s.  They  included  additional 
tomography  results  from  periods  of  65  to  120s.  The  tomography  models  provided  average 
Rayleigh  and  Love  wave  dispersion  curves  for  each  2°  x  2°  block  in  the  region.  They  also  used 
the  results  to  determine  the  velocity  structure  by  fitting  the  synthetic  curves  from  simplified  crust 
and  upper  mantle  models  to  the  tomographic  data  for  each  block  via  a  grid  search.  Finally,  they 
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compared  and  contrasted  their  preferred  velocity  model  with  a  number  of  other  published  models 
for  the  region  and  found  that  they  were  similar. 

Brown  et  al.  (2004)  presented  a  new  3D  density  model  for  the  Libya  region  that  included 
six  individual  subsurface  layers  down  to  a  mantle  depth  of  50  km.  Their  model  was  derived  by 
utilizing  over  6000  gravity  measurements,  in  conjunction  with  other  geophysical  and  geological 
data  obtained  from  oil  exploration,  published  research,  and  other  academic  institutions  to 
develop  a  3-D  density  model  for  the  Libya  region.  They  evaluated  their  model  by  using  it  to 
compute  predicted  Free  Air  gravity  valued  for  the  region  and  comparing  their  results  to  the 
observed  Free  Air  gravity  anomalies.  Their  results  are  shown  in  Figure  21 . 

Geophysical  Model  Development 

Our  3D  velocity  model  (LIBYA3D)  consists  of  five  individual  lithospheric  layers.  The 
layers  are  the  uppermost  layer  (a  combination  of  the  Cenozoic  and  the  Mesozoic  layers),  the 
Paleozoic  layer,  the  Precambrian  layer,  the  lower  crust,  and  the  upper  mantle.  Data  for  defining 
each  layer  were  developed  using  previous  3D  geophysical  studies  as  a  priori  information.  We 
relied  heavily  on  a  density  model  developed  by  Brown  et  al  (2004)  as  an  initial  starting  model. 
The  topography  of  the  surface  of  Libya  was  extracted  from  GTOPO30,  which  was  released  by 
the  USGS  (EROS  data  collection,  1996).  The  surficial  sediments  were  combined  with  other 
rocks  of  Cenozoic  and  Mesozoic  age  for  the  first  layer  of  the  model  since  there  was  no  easy  way 
of  separating  these  layers.  The  Nafe-Drake  relationship  was  used  to  assign  a  Vp  of  5.45  km/s  and 
a  Vs  of  3.1  km/s  to  the  uppermost  layer  (Figure  24  &  Figure  25)  of  our  3D  velocity  model  as 
derived  from  density.  Our  assigned  Vp  and  Vs  are  comparable  to  those  of  previous  studies 
conducted  by:  Dial  (1998)  Vp  =  4.75  km/s,  Vs  =  2.6  km/s;  Gumper  and  Pomeroy  (1970)  Vp  =  5.8 
km/s.  Vs  =  3.4  km/s;  Qiu  et  al.  (1996);  Vp  =  5.0  km/s,  Vs  =  2.8  km/s  and  Blanco  and  Spakman 
(1991)  Vp  =  6.2  km/s.  The  IASPEI91  (Kennett,  and  Engdahl,  1991)  model  is  a  global  model  and 
therefore  the  Vp  and  Vs  for  the  shallowest  layers  were  combined  and  remained  undifferentiated. 
The  bottom  of  LIBYA3D’s  upper  layer  was  derived  from  a  1:  2,000,000  structure  contour  map 
of  the  Libyan  Arab  Republic  and  adjacent  areas  (Goudarzi,  and  Smith,  1978).  The  contour  map 
was  drawn  at  intervals  of  250  m.  The  average  thickness  of  the  combined  uppermost  layer  is 
approximately  1  km  with  a  maximum  thickness  of  4.9  km  occurring  in  the  Sirte  Basin  area. 
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Figure  24.  Average  ID  velocity  model  beneath  Libya.  The  wavy  lines  indicate  surfaces  that 
were  digitized  from  data.  A  uniform  depth  of  19  km  was  assumed  for  the  top  of  the  lower 
crust.  Density  values  are  given  in  g/cc. 
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Figure  25.  The  velocity/density  model  resulting  from  our  study.  Dashed  lines  represent  the 
IASPEI91  model. 

The  base  of  the  Mesozoic  layer  is  the  top  of  the  Paleozoic  layer.  Paleozoic  sedimentary 
rocks  outcropping  in  Libya  also  serve  as  good  control  to  determine  the  extent  of  the  Paleozoic 
layer.  We  used  the  Nafe-Drake  relationship  to  assign  a  Vp  of  5.53  km/s  and  a  Vs  of  3.2  km/s  to 
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the  Paleozoic  Layer.  These  velocities  are  again  comparable  to  those  of  previous  studies  from 
Dial  (1998),  Gumper  and  Pomeroy  (1970),  Qiu  et  al.  (1996)  and  Blanko  and  Spakman  (1991). 
The  Paleozoic  layer  has  a  maximum  thickness  of  4.2  km  beneath  the  Hamra  basin,  northwest 
Libya.  The  average  thickness  of  the  Paleozoic  layer  is  1 .06  km.  The  Paleozoic  layer  is  intruded 
by  volcanic  igneous  rocks  in  sections  of  the  Fezzan  Uplift  region,  western  Libya. 

The  top  of  the  Precambrian  basement  rock  was  also  derived  from  the  1:  2,000,000 
structure  contour  map  of  the  Libya  region  (Goudarzi,  and  Smith,  1978)  that  was  drawn  at 
intervals  of  250  m.  The  Vp  and  Vs  for  this  layer  were  set  at  6.2  km/s  and  3.6  km/s,  respectively, 
from  the  Nafe-Drake  relationship.  These  estimates  are  reasonable  extrapolations  of  the  velocity 
for  rocks  of  Precambrian  age  (Dial  1998).  The  average  thickness  of  the  Precambrian  layer  for  the 
study  area  beneath  Libya  is  17.4  km,  with  a  region  of  maximum  thickness  exceeding  19  km  to 
the  west  of  the  Fezzan  Uplift  where  it  outcrops. 

Not  much  information  is  available  for  the  boundary  between  the  Precambrian  basement 
and  the  lower  crust  and  for  this  reason  we  assumed  a  flat  surface  separating  these  two  layers.  We 
used  results  from  previous  studies  (Dial,  1998)  to  set  this  boundary  at  19  km  and  a  density  of  3 
g/cc.  For  this  density  Vp  and  Vs  extrapolated  from  the  Nafe-Drake  relationship  were  7.16  km/s 
and  4.04  km/s,  respectively.  The  thickness  of  the  lower  crust  ranged  from  20.5  km  in  southern 
Libya,  thinning  to  8.4  km  in  the  north. 

The  Moho  surface  reflects  the  thinning  of  the  crust  towards  the  north  of  Libya  as  shown 
by  Dial  (1998)  and  Yousef  (1986).  The  maximum  depth  chosen  for  the  original  density  model 
was  50  km,  as  it  was  greater  than  the  thickest  region  of  the  crust.  Vp  and  Vs  for  this  portion  of 
the  upper  mantle  were  set  at  7.85  km/s  and  4.54  km/s,  respectively  (Figure  24  &  Figure  25). 

P-wave  velocities  of  the  LIBYA3D  velocity  model  at  depths  of  0  km,  1  km,  and  3  km  are 
shown  in  Figure  26.  These  slices  demonstrate  lateral  velocity  variations  within  the  model. 
Specific  features  in  the  velocity  model  can  be  correlated  with  specific  geologic  or  tectonic 
elements  in  Libya.  For  instance,  a  low  velocity  region  indicating  sediment  cover  in  the  northern 
part  of  the  country  is  clearly  visible  in  Figure  26a  and  Figure  26b.  This  low  velocity  region  is 
due  to  the  presence  of  the  Sirte  Basin  in  the  north  and  north  east  of  the  country  and  the  Hamra 
and  Ghadames  basin  to  the  northwest  of  the  country  (Brown  et  al.  2004).  The  oval  shaped  low 
velocity  zone  in  the  south  east  part  of  the  country  corresponds  to  the  Murzuk  basin  (Brown  et.  al, 
2004).  Two  patches  of  high  velocity  in  Figure  26a  correlate  with  the  Fezzan  uplift  in  the  north 
and  the  Nubian  uplift  to  the  south  (Brown  et  al,  2004). 

Overall  the  resulting  velocity/density  model  showed  a  good  correlation  with  the 
IASPEI91  model  in  all  three  categories  (Vp,  Vs  and  density).  The  values  for  the  P-wave  5-  wave 
and  velocities  from  the  IASPEI91  model  slightly  overestimate  velocities  derived  from  densities 
for  subsurface  Libya.  For  the  Moho  both  models  returned  identical  values  for  P-  and  5-wave 
velocities;  however,  a  density  of  3.3  g/cc  was  assigned  to  our  model  as  opposed  to  3.396  g/cc  in 
the  IASPEI91  model.  The  widest  margin  of  disparity  occurred  within  the  lower  crust  where  the 
1ASPEI  91  P-wave  velocity  estimate  is  6.5  km/s.  Other  studies  (Seber  et  al  1996,  McNamara  et 
al.  1996)  also  suggest  that  a  higher  P- wave  velocity  exists  in  North  Africa.  Dial’s  study  led  him 
to  conclude  that  a  value  of  7.16  km/s  is  a  more  reliable  P- wave  velocity  for  the  lower  crust  in 
this  region.  The  IASPEI91  model  shows  no  differentiation  of  the  layers  above  the  Precambrian 
Basement. 
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Figure  26.  Slices  through  the  LIBYA3D  velocity  model  at  (a)  0  km  depth,  (b)  1  km  depth, 
and  (c)  3  km  depth.  Blue  represents  areas  with  a  velocity  (Vp)  of  5.45  km/s  and  (Vs)  of  3.1 
km/s.  Green  represents  areas  of  velocity  (Vp)  of  5.53  km/s  and  (Vs)  of  3.2km/s.  Red 
represents  areas  of  velocity  (Vp)  of  6.2km/s  and  (Vs)  of  3.6  km/s. 


For  our  LIBYA3D  model  however,  having  digitized  the  surface  of  the  top  of  the 
Precambrian  basement  rock  and  the  base  of  the  Mesozoic  surface,  we  were  able  to  assign 
boundaries  to  the  uppermost  layers  of  the  model.  In  the  Libya  region,  the  crust  thins  towards  the 
north  and  so  there  is  significant  varying  thickness  of  the  Moho.  Figure  24  is  an  average 
representation  of  the  1 D  model  so  we  used  the  average  thickness  of  the  upper  mantle  which  is  35 
km,  the  same  as  the  value  used  in  the  IASPEI91  model.  In  the  next  section  will  verify  if  our 
velocity  model  provides  significant  improvement  over  IASPEI91  in  waveform  modeling. 
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ID  Validation  of  3D  Velocity  Model 

We  carried  out  preliminary  seismic  evaluations  on  the  average  1 D  velocity  model  using 
the  velocities  assigned  to  the  individual  layers  as  shown  in  Figure  24.  For  this  evaluation  we 
chose  events  921012  (YYMMDD)  (5.8mw)  and  event  9980528  (5.5mw)  (Figure  22).  We  obtained 
locations  for  these  events  using  a  bootstrap  technique  after  Petroy  and  Weins  (1989).  These 
events  had  large  enough  magnitudes  to  allow  good  propagation  of  seismic  energy  across  the 
Libyan  lithosphere.  Special  care  had  to  be  taken  when  choosing  source-receiver  pathways  since 
the  border  between  the  African  plate  and  the  Eurasian  plate  is  situated  to  the  north  of  Libya.  In 
this  area  there  is  considerable  thinning  of  the  crust.  We  utilized  the  TauP  Toolkit  (after  Crotwell, 
2000)  to  compute  a  predicted  travel  time  for  the  P-wave  and  S-wave  of  the  two  events  at  the 
Tamanrasset  (TAM)  station,  Algeria  and  compared  the  results  to  the  actual  recorded  waveforms. 
The  results  are  shown  in  Figure  27.  The  red  line  shows  the  predicted  P-wave  arrival  time  (To)  for 
the  event  generated  from  our  model.  The  predicted  arrival  times  were  accurate  to  within  1.8 
seconds  for  the  P-wave  arrival  times  and  to  within  5  seconds  in  predicting  the  S-wave  arrival 
times,  proving  that  our  velocity  model  has  validity  in  the  1 D  case. 

We  then  calculated  individual  path  specific  velocity  models  for  the  same  two  events  in 
order  to  analyze  our  1.5D  velocity  model.  This  was  done  by  calculating  the  average  thickness  of 
each  layer  along  a  straight  path  from  the  event  source  to  TAM  and  assigning  the  resulting 
thickness  to  that  layer.  This  is  different  from  the  average  1 D  layer  thickness,  as  the  model  is  path 
specific.  The  synthetic  waveform  was  calculated  based  on  the  reflectivity  technique  presented  by 
Kennett  (1983)  as  implemented  by  Randall  (1985).  We  used  the  Seismic  Analysis  Code  (SAC) 
(Goldstein,  1 999)  to  carry  out  the  seismic  analysis  for  this  portion  of  the  study.  The  results  of  the 
modeling  of  the  two  events  are  shown  in  Figure  28.  The  seismograms  are  filtered  using  a  band 
pass  Butterworth  filter  of  10- 100s  as  the  most  significant  waveform  information  is  found  in  this 
frequency  range.  For  event  921012  we  observed  a  close  correlation  between  the  synthetic  and  the 
recorded  seismogram  on  the  vertical  and  transverse  components,  while  we  did  not  receive  a  good 
fit  on  the  radial  component.  For  event  980528  we  again  observed  favorable  correlation  on  the 
vertical  component  but  mismatched  the  radial  and  transverse  components  (Figure  28). 
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Figure  27.  Unfiltered  seismogram  of  a)  event  921012  and  b)  event  980528  recorded  at 
TAM.  The  arrival  time  of  the  P-wave  as  predicted  by  our  ID  model  (red  lines)  and  as  well 
as  the  predicted  S-wave  arrival  time  (blue  lines). 
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Error  Analysis 

We  compared  our  predicted  waveforms  to  the  data  by  calculating  the  rms  of  the 
mismatches  between  the  observed  and  predicted  waveforms  for  events  921012  and  980528 
recorded  at  TAM.  These  events  were  used  as  the  waveforms  had  suitably  large  enough 
amplitudes  to  be  modeled.  The  waveforms  were  calculated  as  a  function  of  depth  to  validate  the 
initial  depth  provided  by  the  Harvard  CMT.  The  initial  depths  were  found  to  be  the  most  suitable 
depth  for  performing  the  waveform  modeling.  We  then  filtered  the  waveforms  at  different  pass 
bands  and  the  result  of  comparisons  done  at  the  10- 100s  pass  band  is  shown  in  Figure  29  and 
Table  1.  The  comparison  was  performed  over  the  first  10  seconds  of  the  wave  train  and  Figure 
29  shows  the  resulting  error  analysis  for  the  three  different  components  of  the  two  sample  events 
using  the  L1BYA3D,  the  global  IASPEI91  model  and  the  15  layer  model  proposed  by  Dial 
(1998). 

In  all  cases  the  LIBYA3D  model  returned  a  more  accurately  predicted  waveform  for  the 
sampled  seismic  events  than  the  IASPEI91  model.  This  was  not  the  case  when  LIBYA3D  was 
compared  to  the  Dial  15  layer  model.  Dial’s  (1998)  model  was  derived  from  iterative  seismic 
waveform  modeling  of  predicted  waveforms  compared  against  earthquake  data.  This  analysis  of 
the  waveform  data  suggests  that  Dial’s  (1998)  model  may  have  better  accounted  for  surficial 
layers.  However,  until  we  are  able  to  perform  2D  and  3D  modeling  we  may  not  be  able  to  see  the 
true  improvement  of  LIBYA3D  as  compared  to  the  other  models. 
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Figure  28.  Waveform  comparison  between  synthetic  and  data  for  events  921012  and 
980528  used  in  the  validation  study.  The  seismograms  were  recorded  at  TAM.  A  band  pass 
(10  -  100s)  Butterworth  filter  has  been  applied  to  the  data  and  the  synthetic  seismograms. 
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Figure  29.  Plots  of  the  results  of  our  1.5D  error  analysis  exercise.  The  seismograms  were 
recorded  at  TAM. 
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Table  4.  Comparison  of  error  analysis  results  (in  amplitudes  of  meters)  of  LIBYA3D, 
1ASPE191  and  the  Dial  Model. 


Component 

Velocity  Model 

921020 

libya3D 

IASPEI91 

Dial (1998) 

vertical 

radial 

transverse 

6.20E-04 

3.25E-05 

2.25E-06 

8.25E-04 

5.55E-05 

2.40E-06 

5.25E-04 

1.70E-05 

2.10E-06 

980528 

libya3D 

IASPEI91 

Dial (1998) 

vertical 

radial 

transverse 

3.45E-04 

3.26E-05 

1.39E-04 

4.25E-04 

3.73E-05 

2.30E-04 

2.85E-04 

3.71E-05 

2.30E-04 

Conclusions 

This  paper  presents  a  new  3D  velocity  model  (LIBYA3D)  for  the  Libya  region. 
LIBYA3D  was  derived  from  a  3D  density  model  after  Brown  et  al  (2004)  and  is  similar  in  many 
ways  to  other  previous  regional  and  global  velocity  models.  Tests  results  obtained  by  using  an 
average  ID  model  derived  from  the  3D  structure  validated  our  model  as  being  effective  in 
predicting  P-  and  S-wave  arrival  times.  We  were  able  to  further  validate  our  model  by  producing 
regional  synthetic  waveforms  that  resembled  the  actual  waveforms.  In  error  analysis  tests  of 
LIBYA3D  using  1.5D  approximations  to  the  path  returned  error  values  that  were  comparable  to 
other  regional  velocity  models  developed  from  actual  seismic  waveform  modeling,  while 
outperforming  the  global  IASPEI91  model. 

Due  to  lack  of  events  occurring  in  the  region  with  suitable  source  receiver  paths,  we  were 
unable  to  thoroughly  test  our  model.  However,  we  remain  convinced  that  our  3D  velocity  model 
is  unique  for  the  region  and  it  is  potentially  an  excellent  starting  model  for  anyone  interested  in 
carrying  out  inversion  involving  a  velocity  model  for  the  region.  Developing  such  a  3D  model 
also  improves  our  ability  for  CTBT  monitoring  in  a  complex  geologic  region  such  as  Libya. 
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Section  5 


ANALYSIS  OF  HIGH  FREQUENCY  PROPAGATION  IN  THE  LIBYAN  REGION 
UTILIZING  A  3-D  GEOPHYSICAL  MODEL 

Summary 

We  utilize  a  3-D  geophysical  model,  “LIBYA3D,”  in  the  context  of  regional 
discrimination  to  investigate  small  to  moderate  size  natural  events  (m  3.7  -  5.7)  in  the  Libyan 

b 

region  occurring  from  May  1990  to  May  2000.  Observations  of  Pn  and  Sn  energy  at  varying 
frequency  ranges  (0.5  -  8  Hz)  and  propagation  paths  recorded  by  broadband  instrumentation  led 
us  to  investigate  the  attenuation  effects  of  the  geologic  structure  of  the  Libyan  region  and  the 
adjacent  Mediterranean  Sea.  We  perform  multi-frequency  analysis  to  discern  which  frequencies 
have  the  optimum  signal  to  noise  ratio  (SNR)  for  these  regional  phases.  Pn  and  Sn  SNR  show 
stronger  amplitudes  at  higher  frequencies  (4-8  Hz)  along  propagation  paths  where  continental 
crust  is  limited.  Conversely,  we  only  observe  strong  Pn  and  Sn  amplitudes  at  lower  frequencies 
(0.5  -  4  Hz)  along  propagation  paths  where  continental  crust  is  prevalent.  We  utilize  LIBYA3D 
to  develop  relationships  in  our  data  as  a  function  of  crustal  path,  Moho  variance,  and  volume. 
These  empirical  observations  could  be  used  to  aid  future  clandestine  explosion  discrimination 
work  in  the  Libyan  region  if  the  need  arises. 

Introduction 

To  monitor  a  zero  yield  nuclear  explosion  treaty,  we  must  be  able  to  analyze  events  of  all 
magnitudes  that  occur  throughout  the  world.  However,  it  is  difficult  to  analyze  the  effectiveness 
of  discriminating  between  earthquakes  and  explosions  when  there  are  a  small  amount  of  events 
to  analyze.  Compounding  problems  arise  with  small  magnitude  events  recorded  at  regional 
distances,  since  path  effects  can  lead  to  complex  regional  seismograms  causing  misclassification 
of  events  (e.g.,  Velasco  et  al.,  2003).  Many  regions  in  the  world  exhibit  these  attributes  of  few 
events  and  complex  regional  propagation,  forcing  new  approaches  to  be  developed  in  order  to 
effectively  monitor  for  underground  nuclear  tests. 

The  Libyan  region  is  a  stable  continental  margin  and  exhibits  some  of  these  attributes, 
where  there  are  few  events,  and  those  that  occur  are  generally  small  in  magnitude.  We  focus  on 
studying  the  effectiveness  of  regional  discrimination  techniques  in  the  Libya  region.  We  utilize 
seismic  events  in  and  around  Libya  (Figure  30,  Table  5),  yet  there  are  no  explosion  data  for 
Libya.  We  are  further  limited  by  the  small  number  of  reliable  stations  at  regional  distances. 
Although  there  are  regional  networks  in  place  in  Egypt  (Helwan)  and  Algeria  (Alger- 
Bouzareah),  the  data  are  not  readily  available  to  us.  To  overcome  these  issues,  we  utilize  a  3-D 
geophysical  model  of  the  Libyan  region,  known  as  “LIBYA3D”  (Brown  et  al.,  2004),  to  help 
improve  regional  discrimination  in  the  region.  We  use  LIBYA3D  to  investigate  path  effects 
associated  with  Signal  to  Noise  Ratio  (SNR).  We  also  identify  path  dependent  cross  spectral 
ratios  (e.g., Walter  et  al.,  1995)  which  show  the  elevated  SNR  frequency  ranges  at  specific 
stations  as  solid  parameters  for  future  verification  work. 
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Figure  30.  Map  of  study  areas  show  earthquakes  as  black  diamonds,  stations  as  black 
triangles,  and  red  great  circle  arcs  represent  source-receiver  paths. 

LIBYA3D  Velocity  Model 

The  LIBYA3D  velocity  model  was  derived  from  over  6000  gravity  measurements,  in 
conjunction  with  other  geophysical  and  geological  data  obtained  from  oil  exploration,  published 
research,  and  other  academic  institutions  to  develop  a  3D  density  model  for  the  Libya  region 
(Brown  et  al.,  2004).  The  3-D  velocity  model,  LIBYA3D  was  derived  from  this  density  model. 
The  model  consists  of  5  lithospheric  layers  reflecting  the  geology  at  depth,  which  includes 
Cenozoic/Mesozoic,  and  Paleozoic  sedimentary  layers  overlying  a  Precambrian  basement.  An 
extra  layer  was  added  to  account  for  the  Mediterranean  Sea.  In  the  Libya  region,  the  crust  thins 
towards  the  north,  and  thus  the  model  features  a  varying  crustal  thickness  as  determined  from 
Dial  (1998)  (Figure  31).  Overall  the  resulting  velocity/density  model  shows  a  good  correlation 
with  the  IASPEI91  model  in  P-  and  S-wave  velocity  and  density  (Kennett  et  al.,  1991).  The 
major  differences  between  IASPEI91  and  LIBYA3D  model  lie  on  a  localized  scale. 


49 


Table  5.  Earthquake  Source  Parameters  (International  Seismic  Centre) 


Date 

Time 

Lat. 

Lon. 

Dep. 

mb 

2000/05/21 

22:45:38.3 

31.1570 

15.7340 

23.6 

3.8 

1990/05/18 

18:27:51.2 

31.7130 

24.7960 

10.0 

4.3 

1990/05/25 

15:44:45.9 

37.8702 

21.0578 

1.0 

3.8 

1990/07/30 

17:52:37.9 

34.4845 

25.5438 

45.8 

4.4 

1991/02/02 

05:34:47.3 

35.1190 

22.1571 

10.0 

3.7 

1991/02/08 

13:50:40.1 

35.2852 

23.0458 

63.8 

4.4 

1991/05/28 

18:58:46.6 

37.0235 

22.4817 

20.8 

4.5 

1991/08/08 

22:46:10.9 

35.6467 

11.6387 

22.0 

4.7 

1992/07/01 

15:13:32.8 

31.4246 

15.8143 

10.0 

4.2 

1992/10/12 

13:09:55.4 

29.7566 

31.1360 

21.5 

5.8 

1993/09/09 

00:32:35.5 

32.2394 

20.3859 

10.0 

4.8 

1995/07/22 

06:24:45.5 

34.7502 

23.3754 

72.2 

4.1 

1995/08/02 

01:42:20.6 

33.4384 

25.3153 

33.0 

3.4 

1995/08/27 

19:42:14.7 

38.2610 

15.2912 

25.2 

3.6 

1995/11/25 

14:10:46.0 

28.8123 

34.9227 

14.0 

3.8 

1996/01/06 

10:38:23.9 

32.9867 

23.6555 

1.0 

1996/03/23 

20:40:46.9 

36.9910 

21.8653 

30.3 

3.5 

1996/05/19 

14:15:14.2 

34.1160 

25.3398 

37.6 

4.0 

1996/06/21 

13:03:38.5 

37.8847 

21.2394 

91.1 

3.5 

1996/10/09 

16:57:40.1 

36.8145 

21.2931 

13.6 

4.0 

1996/10/09 

17:01:49.0 

34.4148 

32.1809 

50.2 

4.1 

1996/10/30 

10:43:02.9 

32.4343 

20.5326 

10.0 

4.6 

1997/02/01 

23:08:36.0 

35.0318 

24.4391 

10.0 

3.9 

1997/05/10 

14:03:03.6 

24.8180 

24.4750 

0.0 

3.9 

1998/01/13 

12:34:21.1 

37.3992 

20.7101 

15.1 

3.6 

1998/04/12 

22:06:31.0 

37.4186 

20.7159 

10.0 

4.4 

1998/05/28 

18:33:28.4 

31.4459 

27.6426 

10.0 

5.5 

1998/08/08 

16:39:28.7 

36.2939 

22.7695 

5.0 

3.5 

1999/03/31 

09:42:16.1 

32.9550 

14.8380 

32.2 

4.3 

1999/06/17 

10:18:02.2 

32.5410 

24.0530 

10.0 

4.3 

1999/11/07 

14:43:36.0 

32.6600 

21.1740 

10.0 

4.4 

Figure  31.  A  portion  of  the  LIBYA3D  Moho  boundary  DEM  used  in  our  study  showing 
crustal  thickness  under  Libya,  note  the  thinning  crust  to  the  north  as  Libya  approaches  the 
Mediterranean  Sea. 


Regional  Seismic  Data 

We  gathered  seismic  data  for  events  that  occurred  from  18  May  1990  to  21  May  2000 
in  the  Libyan  region  and  the  adjacent  Mediterranean  Sea.  Event  times  and  locations  were 
obtained  from  the  International  Seismic  Centre  (ISC).  We  have  data  consisting  of  3 1  natural 
events,  7  located  in  Libya,  1  located  in  Egypt,  with  the  remainder  located  in  the  adjacent  area  of 
the  Mediterranean  Sea.  Magnitudes  range  from  (m  )  3.7  to  5.7. 

b 

Waveform  data  were  acquired  from  the  Incorporated  Research  Institutions  for 
Seismology’s  (IRIS)  Data  Management  Center  (DMC).  The  stations  used  in  this  research  effort 
are  Kottamia,  Egypt  (KEG),  Wield  Dalam,  Malta  (WDD),  Tamanrasset,  Algeria  (TAM),  Anoyia, 
Crete  (IDI),  and  Gafsa,  Tunisia  (GFA),  all  with  broadband  instrumentation.  We  remove 
instrument  response  and  set  origin  times  for  each  event.  Waveform  data  are  then  cut  and  rotated 
to  their  radial  and  transverse  components,  and  finally,  we  apply  a  Butterworth  filter  to  view  the 
waveforms  at  several  frequency  ranges  (0.5-1,  1-2,  2-4,  4-6,  and  6-8  Hz)  (Figure  32).  We  then 
pick  regional  phases  Pn  and  Sn  manually  and  mark  the  header  for  each  waveform. 

We  show  the  02  February  1991  (910202)  event  (Figure  32)  recorded  at  KEG  (A  =  9.67  ; 

O 

AZ  =  120  )  (see  Figure  1);  note  the  dramatic  increase  in  SNR  at  the  higher  frequencies.  This  led 
us  to  investigate  attenuation  effects  of  source-receiver  paths  as  related  to  their  orientation  across 
the  varying  topography  and  bathymetry  of  Libya  and  the  Mediterranean  Sea.  Plotting  all 
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seismograms  at  KEG  and  TAM  for  regional  events  shows  consistent  high  frequency  propagation 
of  Pn  and  Sn  energy  at  KEG  and  low  frequency  propagation  of  Pn  and  Sn  energy  at  TAM. 
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Figure  32.  Vertical  component  multi-frequency  bandpassed  seismogram  of  the  02 
February  1991  event  recorded  at  KEG.  Note  that  SNR  increases  at  the  higher  frequency 
ranges.  This  is  a  good  example  that  led  us  to  question  what  path  effects  were  responsible 
for  the  observed  high  frequency  propagation. 


Two  events  of  particular  interest,  due  to  the  fact  they  were  recorded  at  two  receivers,  are 

the  31  March  1999  (990331)  event  recorded  at  IDI  (A  =  8.66  ;  AZ  =  71.9  )  and  TAM  (A  =  13.1  ; 

AZ  =  221°),  and  the  01  July  1992  (920701)  event  recorded  at  KEG  (A  =  13.9°;  AZ  =  92.0°)  and 

TAM  (A  =  12.6  ;  AZ  =  229  )  (Figure  33)  (see  Figure  30).  The  SNR  for  both  events  at  TAM  is 
low,  from  4  to  8  Hz;  however  at  IDI  and  KEG  the  SNR  is  still  larger,  if  not  dominant,  from  4  to 
8  Hz.  Note  the  receiver  distance  to  KEG  is  greater  than  the  receiver  distance  to  TAM  for  920701 
(Figure  33B),  but  the  SNR  is  still  strong  enough  to  be  seen  at  KEG  in  the  high  frequency  ranges. 
This  demonstrates  that  receiver  distance  is  not  the  primary  controlling  parameter  for  signal 
attenuation. 

We  calculate  phase  amplitude  and  signal  to  noise  ratio  (SNR)  by  measuring  the 
maximum  displacement  in  a  ten  second  pre-phase  and  phase  window.  By  comparing  the  SNR  of 
Pn  against  frequency  at  a  specific  station  (Figure  34A-C)  we  observe  trends  in  station-frequency 
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dependence.  These  trends  were  normalized  by  dividing  each  point  by  the  SNR  at  the  lowest 
frequency  range  for  each  event  (0.5-1. 0  Hz).  This  highlights  the  difference  in  high  frequency 
propagation.  As  seen  in  these  diagrams,  KEG  conveys  the  most  dramatic  results  with  the  SNR 
peaking  at  6-8  Hz  in  10  of  the  14  events  recorded  there.  IDI  and  WDD  show  a  trend  similar  to 
KEG  with  good  SNR  at  4  to  8  Hz,  but  they  actually  peak  at  3  Hz  for  most  of  these  Intra- 
Mediterranean  events.  We  combine  their  results  due  to  their  similar  trends,  as  well  as  the  fact 
that  their  receiver  environments  are  the  same  topographically.  TAM  shows  expected  high 
frequency  attenuation  of  Pn  energy. 
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Figure  33.  Vertical  component  multi-frequency  seismograms  from  the  990331  and  920701 
events.  Notice  the  different  high-frequency  propagation  when  comparing  IDI  and  KEG  to 
TAM.  Also  note  that  even  though  the  receiver  distance  is  greater  to  KEG  than  to  TAM  for 
920701  a  distinguishable  signal  is  still  present  from  4  to  8  Hz  at  KEG. 
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Figure  34.  SNR-Frequency  diagrams  show  trends  in  station-frequency  dependence,  note 
the  unusual  high  frequency  propagation  displayed  at  KEG  accompanied  by  the  expected 
high  frequency  attenuation  at  TAM. 


Application  of  LIBYA3D 

Using  LIBYA3D  (see  Figure  31)  we  were  able  to  calculate  variable  topographic  statistics 
concerning  the  crust  and  the  Moho  related  to  the  source-receiver  paths.  With  the  use  of  a  0.5 
degree  wide  swath,  we  calculate  the  crustal  volume  for  each  source-receiver  path.  We  also  use  a 
simple,  but  accurate  map  measurement  to  determine  the  percentage  of  continental  crust  in  each 
source-receiver  path.  Finally,  using  LIBYA3D  we  calculate  the  standard  deviation  of  the  Moho, 
which  is  similar  to  Zhang  and  Lay’s  (1994)  surface  roughness  calculation,  of  each  source- 
receiver  path.  We  then  calculate  the  cross-spectral  amplitude  and  cross-spectral  SNR  values  of 
Pn  ( 1-2  /  6-8  Hz)  for  each  event;  in  this  fashion,  we  can  define  the  events  with  high-frequency 
propagation  by  their  low  cross-spectral  value.  We  plot  these  topographic  statistics  against  the 
cross-spectral  values  and  seek  correlations  (Figure  35A-F)  (e.g.,  Zhang  &  Lay,  1994).  Due  to 
the  small  amount  of  data  that  we  have  to  work  with  for  each  plot  we  define  a  correlation 
coefficient  (r)  greater  than  0.5  as  a  strong  correlation. 

At  KEG,  we  find  correlations  between  high  frequency  propagation  and  the  percentage  of 
continental  crust  in  the  source-receiver  path  when  compared  to  cross-spectral  amplitude  and 
cross-spectral  SNR  (Figure  35A,  B)  and  a  weak  correlation  between  high  frequency  propagation 
and  crustal  volume.  For  IDI,  WDD,  and  TAM  there  are  strong  correlations  when  comparing 
percentage  of  continental  crust  and  crustal  volume  to  cross-spectral  SNR  (Figure  35C,  D). 
However,  these  correlations  break  down  when  compared  to  cross-spectral  amplitude.  Finally, 
we  find  two  distinct  correlations  at  KEG  for  standard  deviation  of  the  Moho  (Figure  35E,  F). 
When  compared  to  cross-spectral  SNR  a  positive  slope  is  exhibited,  however,  when  compared  to 
cross-spectral  amplitude  a  negative  slope  is  produced. 
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Figure  35.  Cross-spectral  value  plots  show  a  good  correlation  between  the  continental 
crust  percentage  in  the  source-receiver  path  and  cross-spectral  values  at  KEG  (A  &  B). 
Strong  correlations  are  also  shown  for  IDI,  TAM,  and  WDD  when  compared  to  crustal 
volume  and  continental  crust  percentage  (C  &  D).  Finally,  we  observe  converse  trends  at 
KEG  when  compared  to  the  standard  deviation  of  the  Moho  depth  (E  &  F). 
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Discussion 

A  3-D  velocity  model  can  aid  in  improving  discrimination  between  earthquakes  and 
explosions  by  providing  a  means  for  determining  the  contributing  factors  in  regional  wave 
propagation.  Formal  approaches  to  regional  waveform  propagation,  such  as  Magnitude  Distance 
and  Amplitude  Corrections  (MDAC)  (Taylor  et  al.,  2002)  rely  on  developing  magnitude  and 
amplitude  corrections  based  on  empirical  data  at  various  frequencies.  However,  when  few 
events  exist  due  to  the  nature  of  the  tectonic  region,  these  empirical  based  approaches  break 
down.  Utilizing  a  3D  model  developed  using  other  geophysical  data  may  be  the  only  means  to 
correct  for  high  frequency  regional  propagation.  Although  we  do  not  formally  develop 
corrections,  we  can  extend  this  application  of  the  model  to  develop  crustal  path  corrections  to 
high  frequency  propagation. 

We  see  through  our  various  geophysical  plots  that  there  are  several  parameters  that  affect 
amplitude  and  SNR  at  varying  frequencies  for  a  particular  phase,  in  our  case,  Pn.  By  having  all 
of  our  plots  dependent  on  the  same  cross-spectral  values  we  can  empirically  gauge  which 
parameters  have  the  greatest  effect  on  these  values.  We  identify  a  strong  correlation  between 
percentage  of  continental  crust  and  cross-spectral  SNR  at  KEG.  More  importantly,  the 
correlations  between  crustal  volume  and  Moho  standard  deviation  with  cross  spectral  values  not 
only  show  solid  parameters  for  future  correction  and  verification  work,  but  also  serve  to  validate 
the  LIBYA3D  model. 

We  recognize  that  these  correlations  are  not  perfect,  but  nature  rarely  exhibits  such 
perfect  fits.  This  leads  us  to  further  examine  outlying  events  in  our  plots.  We  see  that  high 
frequency  propagation  is  most  stringently  controlled  by  the  percentage  of  continental  crust  in  the 
source-receiver  path,  which  relates  to  the  crustal  volume.  Secondarily,  standard  deviation 
(roughness)  of  the  Moho  boundary  influences  high  frequency  propagation.  We  also  see  that 
other  source  parameters  such  as  depth,  magnitude,  and  finally  distance  all  play  a  hand  in 
controlling  the  high-frequency  propagation  of  Pn.  Since  we  know  that  velocity  gradients  above 
and  below  the  Moho  trap  seismic  energy  resulting  in  larger  Pn  amplitudes  (Stein  &  Wysession, 
2003)  and  it  has  been  shown  that  Pn  attenuates  less  rapidly  at  higher  frequencies  if  the  spherical 
earth  velocity  gradient  is  zero  or  greater  (Sereno  &  Givens,  1990);  we  feel  that  this  is  also  a 
parameter  controlling  high  frequency  propagation  in  our  events. 

Conclusions 

By  performing  a  multi-frequency  analysis,  we  find  that  Pn  and  Sn  show  stronger 
amplitudes  at  higher  frequencies  (4-8  Hz)  along  propagation  paths  where  continental  crust  is 
limited.  Conversely,  we  only  observe  strong  Pn  and  Sn  amplitudes  at  lower  frequencies  (0.5  -  4 
Hz)  along  propagation  paths  where  continental  crust  is  prevalent.  We  discover  that  at  regional 
station  KEG,  we  see  a  strong  correlation  between  percentage  of  continental  crust  and  signal 
propagation  in  the  high  frequency  ranges,  and  this  relationship  persists  as  well  for  IDI,  WDD, 
and  TAM,  even  though  there  is  a  lack  of  data.  Most  of  the  paths  for  IDI  and  WDD  are  purely 
oceanic,  since  they  are  island  stations,  while  the  paths  for  TAM  are  purely  continental,  which 
skews  the  results.  We  feel  that  given  a  full  range  of  data  this  parameter  will  consistently  predict 
high  frequency  propagation  events  that  would  aid  in  verification  work  in  regions  of  mixed  crusts. 
Using  LIBYA3D,  we  found  correlations  for  all  four  stations  when  comparing  crustal  volumes  to 
cross-spectral  values.  We  feel  more  data  are  needed  before  drawing  any  solid  conclusions  on  the 
relationship  between  crustal  volume  and  high  frequency  propagation. 
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We  find  interesting  results  at  KEG  when  comparing  cross-spectral  values  against  the 
standard  deviation  (roughness)  of  the  Moho.  We  see  that  greater  roughness  degrades  the  SNR  at 
higher  frequencies  (4-8  Hz),  but  simultaneously  causes  less  rapid  attenuation  at  these  higher 
frequencies.  We  question  if  there  is  a  link  between  roughness  of  the  Moho  and  its  velocity 
gradient,  which  is  known  to  increase  the  amplitude  of  the  Pn  phase.  Finally,  we  conclude  that  3- 
D  geophysical  models  such  as  LIBYA3D  are  valuable  for  conducting  regional  amplitude 
correction  and  verification  work  in  regions  where  little  or  no  explosion  data  are  available  This 
becomes  even  more  important  since  this  is  the  case  for  many  stations  currently  monitoring  the 
Comprehensive  Test  Ban  Treaty  (CTBT). 
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